Cod diet model
Max Lindmark & Michele Casini
27 October, 2021
Intro
Fit density (also referred to as CPUE) model with environmental predictors and use that to calculate weighted mean dissolved oxygen, temperature and depth of Baltic cod
# Load libraries, install if needed
library(tidyverse); theme_set(theme_light(base_size = 10))
library(readxl)
library(tidylog)
library(RCurl)
library(viridis)
library(RColorBrewer)
library(patchwork)
library(janitor)
library(icesDatras)
library(mapdata)
library(patchwork)
library(rgdal)
library(raster)
library(sf)
library(rgeos)
library(chron)
library(lattice)
library(ncdf4)
library(marmap)
library(rnaturalearth)
library(rnaturalearthdata)
library(mapplots)
library(EnvStats)
library(qwraps2)
#remotes::install_github("pbs-assess/sdmTMB")
library(sdmTMB)# To load entire cache in interactive r session, do: qwraps2::lazyload_cache_dir(path = "R/analysis/spatial_trend_models_cache/html")For maps
# Specify map ranges
ymin = 55; ymax = 58; xmin = 14; xmax = 20
map_data <- rnaturalearth::ne_countries(
scale = "medium",
returnclass = "sf", continent = "europe")
# Crop the polygon for plotting and efficiency:
# st_bbox(map_data) # find the rough coordinates
swe_coast <- suppressWarnings(suppressMessages(
st_crop(map_data,
c(xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax))))
# Transform our map into UTM 33 coordinates, which is the equal-area projection we fit in:
utm_zone33 <- 32633
swe_coast_proj <- sf::st_transform(swe_coast, crs = utm_zone33)
ggplot(swe_coast_proj) + geom_sf() 
# Define plotting theme for main plot
theme_plot <- function(base_size = 10, base_family = "") {
theme_light(base_size = 10, base_family = "") +
theme(
axis.text.x = element_text(angle = 90),
axis.text = element_text(size = 8),
legend.text = element_text(size = 8),
legend.title = element_text(size = 8),
legend.position = "bottom",
legend.key.height = unit(0.2, "cm"),
legend.margin = margin(0, 0, 0, 0),
legend.box.margin = margin(-5, -5, -5, -5),
strip.text = element_text(size = 8, colour = 'black', margin = margin()),
strip.background = element_rect(fill = "grey90")
)
}
# Define plotting theme for facet_wrap map with years
theme_facet_map <- function(base_size = 10, base_family = "") {
theme_light(base_size = 10, base_family = "") +
theme(
axis.text.x = element_text(angle = 90),
axis.text = element_text(size = 6),
strip.text = element_text(size = 8, colour = 'black', margin = margin()),
strip.background = element_rect(fill = "grey90")
)
}
# Make default base map plot
plot_map_raster <-
ggplot(swe_coast_proj) +
geom_sf(size = 0.3) +
labs(x = "Longitude", y = "Latitude") +
theme_facet_map(base_size = 14)
# Function to convert from lat lon to UTM
LongLatToUTM <- function(x, y, zone){
xy <- data.frame(ID = 1:length(x), X = x, Y = y)
coordinates(xy) <- c("X", "Y")
proj4string(xy) <- CRS("+proj=longlat +datum=WGS84") ## for example
res <- spTransform(xy, CRS(paste("+proj=utm +zone=",zone," ellps=WGS84",sep='')))
return(as.data.frame(res))
}Read data
d <- readr::read_csv("/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/bentfish/data/for_analysis/stomach/full_stomach_data_21.10.25.csv") %>% dplyr::select(-X1)
#> Warning: Missing column names filled in: 'X1' [1]
#>
#> ── Column specification ────────────────────────────────────────────────────────
#> cols(
#> .default = col_double(),
#> Sample = col_character(),
#> Length.code = col_character(),
#> Prey.sp.code = col_character(),
#> Comment = col_character(),
#> Parasites.in.stomach = col_character(),
#> Coun. = col_character(),
#> Ship = col_logical(),
#> Cruise = col_character(),
#> Predator.code = col_character(),
#> Sex = col_character(),
#> Q.year = col_character(),
#> Date = col_date(format = ""),
#> Index = col_character(),
#> Ices.rect = col_character(),
#> source = col_character(),
#> Number = col_logical(),
#> Sample.type = col_logical(),
#> transect = col_logical(),
#> Depthstep = col_logical(),
#> Gonad.weight.roundfish = col_logical()
#> # ... with 6 more columns
#> )
#> ℹ Use `spec()` for the full column specifications.
#> Warning: 94462 parsing failures.
#> row col expected actual file
#> 3364 Ship 1/0/T/F/TRUE/FALSE SVEA '/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/bentfish/data/for_analysis/stomach/full_stomach_data_21.10.25.csv'
#> 3365 Ship 1/0/T/F/TRUE/FALSE SVEA '/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/bentfish/data/for_analysis/stomach/full_stomach_data_21.10.25.csv'
#> 3366 Ship 1/0/T/F/TRUE/FALSE SVEA '/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/bentfish/data/for_analysis/stomach/full_stomach_data_21.10.25.csv'
#> 3367 Ship 1/0/T/F/TRUE/FALSE SVEA '/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/bentfish/data/for_analysis/stomach/full_stomach_data_21.10.25.csv'
#> 3368 Ship 1/0/T/F/TRUE/FALSE SVEA '/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/bentfish/data/for_analysis/stomach/full_stomach_data_21.10.25.csv'
#> .... .... .................. ...... ................................................................................................................
#> See problems(...) for more details.
# Add UTM coordinates
utm_coords <- LongLatToUTM(d$Long, d$Lat, zone = 33)
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
#> prefer_proj): Discarded datum Unknown based on WGS84 ellipsoid in CRS definition
d$X <- utm_coords$X/1000
d$Y <- utm_coords$Y/1000Explore data
head(data.frame(d))
#> SubDiv Year Month Day N.with.food N.regurgitated N.skeletal N.empty HN Sample
#> 1 28 2017 11 27 1 NA NA NA 52 507
#> 2 28 2017 11 27 1 NA NA NA 52 507
#> 3 28 2017 11 27 1 NA NA NA 52 507
#> 4 28 2017 11 27 1 NA NA NA 52 507
#> 5 28 2017 11 27 1 NA NA NA 52 507
#> 6 28 2017 11 27 1 NA NA NA 52 507
#> Length.code Prey.sp.code Prey.size Prey.weight Prey.nr Stage.digestion
#> 1 <NA> Saduria entomon 20 0.21 1 1
#> 2 <NA> Saduria entomon 16 0.10 1 1
#> 3 <NA> Saduria entomon 20 0.15 1 1
#> 4 <NA> Saduria entomon 20 0.17 1 1
#> 5 <NA> Saduria entomon NA 0.27 2 1
#> 6 <NA> Macoma balthica NA 0.19 1 1
#> Comment Parasites.in.stomach Quarter Coun. Ship Cruise Pred.size.mm
#> 1 <NA> <NA> 4 SWE NA BITS 200
#> 2 <NA> <NA> 4 SWE NA BITS 200
#> 3 <NA> <NA> 4 SWE NA BITS 200
#> 4 <NA> <NA> 4 SWE NA BITS 200
#> 5 <NA> <NA> 4 SWE NA BITS 200
#> 6 <NA> <NA> 4 SWE NA BITS 200
#> Pred.weight.g Predator.gutted.weight Predator.code Sex Maturity Gall.content
#> 1 100 92 FLE M 4 NA
#> 2 100 92 FLE M 4 NA
#> 3 100 92 FLE M 4 NA
#> 4 100 92 FLE M 4 NA
#> 5 100 92 FLE M 4 NA
#> 6 100 92 FLE M 4 NA
#> Q.year Age Date Index Lat Long Depth.catch Ices.rect
#> 1 4_2017 3 2017-11-27 OXBH.2017.52 57.46667 19.23333 67 43G9
#> 2 4_2017 3 2017-11-27 OXBH.2017.52 57.46667 19.23333 67 43G9
#> 3 4_2017 3 2017-11-27 OXBH.2017.52 57.46667 19.23333 67 43G9
#> 4 4_2017 3 2017-11-27 OXBH.2017.52 57.46667 19.23333 67 43G9
#> 5 4_2017 3 2017-11-27 OXBH.2017.52 57.46667 19.23333 67 43G9
#> 6 4_2017 3 2017-11-27 OXBH.2017.52 57.46667 19.23333 67 43G9
#> source Number Sample.type transect Depthstep Gonad.weight.roundfish
#> 1 Max Postdoc NA NA NA NA NA
#> 2 Max Postdoc NA NA NA NA NA
#> 3 Max Postdoc NA NA NA NA NA
#> 4 Max Postdoc NA NA NA NA NA
#> 5 Max Postdoc NA NA NA NA NA
#> 6 Max Postdoc NA NA NA NA NA
#> Liver.weight.roundfish Stomach.content Perc.stomac.content comments Processed
#> 1 NA NA NA NA NA
#> 2 NA NA NA NA NA
#> 3 NA NA NA NA NA
#> 4 NA NA NA NA NA
#> 5 NA NA NA NA NA
#> 6 NA NA NA NA NA
#> Unique.pred.id X Y
#> 1 2017_4_52_507_FLE 753.841 6377.247
#> 2 2017_4_52_507_FLE 753.841 6377.247
#> 3 2017_4_52_507_FLE 753.841 6377.247
#> 4 2017_4_52_507_FLE 753.841 6377.247
#> 5 2017_4_52_507_FLE 753.841 6377.247
#> 6 2017_4_52_507_FLE 753.841 6377.247
sort(colnames(d))
#> [1] "Age" "Comment" "comments"
#> [4] "Coun." "Cruise" "Date"
#> [7] "Day" "Depth.catch" "Depthstep"
#> [10] "Gall.content" "Gonad.weight.roundfish" "HN"
#> [13] "Ices.rect" "Index" "Lat"
#> [16] "Length.code" "Liver.weight.roundfish" "Long"
#> [19] "Maturity" "Month" "N.empty"
#> [22] "N.regurgitated" "N.skeletal" "N.with.food"
#> [25] "Number" "Parasites.in.stomach" "Perc.stomac.content"
#> [28] "Pred.size.mm" "Pred.weight.g" "Predator.code"
#> [31] "Predator.gutted.weight" "Prey.nr" "Prey.size"
#> [34] "Prey.sp.code" "Prey.weight" "Processed"
#> [37] "Q.year" "Quarter" "Sample"
#> [40] "Sample.type" "Sex" "Ship"
#> [43] "source" "Stage.digestion" "Stomach.content"
#> [46] "SubDiv" "transect" "Unique.pred.id"
#> [49] "X" "Y" "Year"
# [1] "Age" "Comment" "comments" "Coun."
# [5] "Cruise" "Date" "Day" "Depth.catch"
# [9] "Depthstep" "Gall.content" "Gonad.weight.roundfish" "HN"
# [13] "Ices.rect" "Index" "Lat" "Length.code"
# [17] "Liver.weight.roundfish" "Long" "Maturity" "Month"
# [21] "N.empty" "N.regurgitated" "N.skeletal" "N.with.food"
# [25] "Number" "Parasites.in.stomach" "Perc.stomac.content" "Pred.size.mm"
# [29] "Pred.weight.g" "Predator.code" "Predator.gutted.weight" "Prey.nr"
# [33] "Prey.size" "Prey.sp.code" "Prey.weight" "Processed"
# [37] "Q.year" "Quarter" "Sample" "Sample.type"
# [41] "Sex" "Ship" "source" "Stage.digestion"
# [45] "Stomach.content" "SubDiv" "transect" "Unique.pred.id"
# [49] "X"
# Plot data in space, color by survey
plot_map_raster +
geom_point(data = d, aes(x = X * 1000, y = Y * 1000, color = Cruise), size = 0.5) +
scale_color_brewer(palette = "Set2") +
facet_wrap(~ Cruise, ncol = 2) +
theme_plot()
# Calculate total weight of prey by predator ID and prey species (i.e., across prey sizes)
d_wide <- d %>%
drop_na(Prey.weight) %>%
group_by(Unique.pred.id, Prey.sp.code) %>%
summarise(tot_biom_per_prey = sum(Prey.weight)) %>%
ungroup() %>%
pivot_wider(names_from = Prey.sp.code, values_from = tot_biom_per_prey) %>%
mutate(across(everything(), ~replace_na(.x, 0))) %>% # Replace NA with 0
janitor::clean_names()
#> drop_na: removed 1,216 rows (6%), 20,515 rows remaining
#> group_by: 2 grouping variables (Unique.pred.id, Prey.sp.code)
#> summarise: now 16,767 rows and 3 columns, one group variable remaining (Unique.pred.id)
#> ungroup: no grouping variables
#> pivot_wider: reorganized (Prey.sp.code, tot_biom_per_prey) into (Sprattus sprattus, Bylgides sarsi, Saduria entomon, Clupea harengus, Pisces, …) [was 16767x3, now 8311x113]
#> mutate: changed 7,085 values (85%) of 'Sprattus sprattus' (7085 fewer NA)
#> changed 7,030 values (85%) of 'Bylgides sarsi' (7030 fewer NA)
#> changed 6,431 values (77%) of 'Saduria entomon' (6431 fewer NA)
#> changed 7,616 values (92%) of 'Clupea harengus' (7616 fewer NA)
#> changed 7,380 values (89%) of 'Pisces' (7380 fewer NA)
#> changed 8,185 values (98%) of 'Gadus morhua' (8185 fewer NA)
#> changed 7,956 values (96%) of 'Cumacea' (7956 fewer NA)
#> changed 7,345 values (88%) of 'Limecola balthica' (7345 fewer NA)
#> changed 6,845 values (82%) of 'Mysis mixta' (6845 fewer NA)
#> changed 8,307 values (>99%) of 'Corophium volutator' (8307 fewer NA)
#> changed 7,639 values (92%) of 'Halicryptus spinulosus' (7639 fewer NA)
#> changed 8,203 values (99%) of 'Stone' (8203 fewer NA)
#> changed 7,662 values (92%) of 'Pontoporeia femorata' (7662 fewer NA)
#> changed 8,285 values (>99%) of 'Enchelyopus cimbrius' (8285 fewer NA)
#> changed 7,946 values (96%) of 'Clupeidae' (7946 fewer NA)
#> changed 7,916 values (95%) of 'Gobiidae' (7916 fewer NA)
#> changed 8,109 values (98%) of 'Neomysis integer' (8109 fewer NA)
#> changed 7,989 values (96%) of 'Monoporeia affinis' (7989 fewer NA)
#> changed 7,865 values (95%) of 'Gammarus sp.' (7865 fewer NA)
#> changed 8,034 values (97%) of 'Mysidae' (8034 fewer NA)
#> changed 8,265 values (99%) of 'Scoloplos armiger' (8265 fewer NA)
#> changed 8,297 values (>99%) of 'plastic' (8297 fewer NA)
#> changed 8,216 values (99%) of 'Priapulus caudatus' (8216 fewer NA)
#> changed 8,257 values (99%) of 'Bivalvia' (8257 fewer NA)
#> changed 8,295 values (>99%) of 'Hyperia galba' (8295 fewer NA)
#> changed 8,091 values (97%) of 'Mytilus sp.' (8091 fewer NA)
#> changed 8,019 values (96%) of 'Crangon crangon' (8019 fewer NA)
#> changed 8,263 values (99%) of 'Idotea balthica' (8263 fewer NA)
#> changed 8,309 values (>99%) of 'Trachurus trachurus' (8309 fewer NA)
#> changed 8,063 values (97%) of 'Gasterosteus aculeatus' (8063 fewer NA)
#> changed 8,309 values (>99%) of 'Annelida' (8309 fewer NA)
#> changed 8,221 values (99%) of 'Algae' (8221 fewer NA)
#> changed 8,296 values (>99%) of 'Platichthys flesus' (8296 fewer NA)
#> changed 8,307 values (>99%) of 'Priapulidae' (8307 fewer NA)
#> changed 8,306 values (>99%) of 'Waste' (8306 fewer NA)
#> changed 8,278 values (>99%) of 'Praunus flexuosus' (8278 fewer NA)
#> changed 8,301 values (>99%) of 'Unidentified mass' (8301 fewer NA)
#> changed 8,308 values (>99%) of 'Merlangius merlangus' (8308 fewer NA)
#> changed 8,252 values (99%) of 'Scales' (8252 fewer NA)
#> changed 8,308 values (>99%) of 'Gadidae' (8308 fewer NA)
#> changed 8,257 values (99%) of 'Crustacea' (8257 fewer NA)
#> changed 8,211 values (99%) of 'Amphipoda' (8211 fewer NA)
#> changed 8,258 values (99%) of 'Spine' (8258 fewer NA)
#> changed 8,309 values (>99%) of 'Pleuronectes platessa' (8309 fewer NA)
#> changed 8,310 values (>99%) of 'Anguilla anguilla' (8310 fewer NA)
#> changed 8,310 values (>99%) of 'Empty' (8310 fewer NA)
#> changed 8,310 values (>99%) of 'Mollusca' (8310 fewer NA)
#> changed 8,309 values (>99%) of 'Pungitius pungitius' (8309 fewer NA)
#> changed 7,462 values (90%) of 'Diastylis rathkei' (7462 fewer NA)
#> changed 8,310 values (>99%) of 'Terebellides stroemii' (8310 fewer NA)
#> changed 8,293 values (>99%) of 'Palaemon sp.' (8293 fewer NA)
#> changed 8,305 values (>99%) of 'Ammodytes tobianus' (8305 fewer NA)
#> changed 7,506 values (90%) of 'NA' (7506 fewer NA)
#> changed 8,309 values (>99%) of 'Cottidae' (8309 fewer NA)
#> changed 8,303 values (>99%) of 'Hediste diversicolor' (8303 fewer NA)
#> changed 8,296 values (>99%) of 'Palaemon elegans' (8296 fewer NA)
#> changed 8,308 values (>99%) of 'Spinachia spinachia' (8308 fewer NA)
#> changed 8,304 values (>99%) of 'Zoarces viviparus' (8304 fewer NA)
#> changed 8,286 values (>99%) of 'Ammodytidae' (8286 fewer NA)
#> changed 8,304 values (>99%) of 'Caridea' (8304 fewer NA)
#> changed 8,310 values (>99%) of 'Amphibalanus improvisus' (8310 fewer NA)
#> changed 8,310 values (>99%) of 'Myoxocephalus quadricornis' (8310 fewer NA)
#> changed 8,305 values (>99%) of 'Phyllodocida' (8305 fewer NA)
#> changed 8,226 values (99%) of 'Remains' (8226 fewer NA)
#> changed 8,302 values (>99%) of 'Palaemonidae' (8302 fewer NA)
#> changed 8,298 values (>99%) of 'Carcinus maenas' (8298 fewer NA)
#> changed 8,310 values (>99%) of 'Gastropoda' (8310 fewer NA)
#> changed 8,300 values (>99%) of 'Sand' (8300 fewer NA)
#> changed 8,308 values (>99%) of 'Cerastoderma glaucum' (8308 fewer NA)
#> changed 8,302 values (>99%) of 'Hyperoplus lanceolatus' (8302 fewer NA)
#> changed 8,285 values (>99%) of 'Mya arenaria' (8285 fewer NA)
#> changed 8,298 values (>99%) of 'Hydrobia sp.' (8298 fewer NA)
#> changed 8,304 values (>99%) of 'Wood' (8304 fewer NA)
#> changed 8,308 values (>99%) of 'Pleuronectiformes' (8308 fewer NA)
#> changed 8,310 values (>99%) of 'Scophthalmus maximus' (8310 fewer NA)
#> changed 8,283 values (>99%) of 'Polychaeta' (8283 fewer NA)
#> changed 8,213 values (99%) of 'Priapulida' (8213 fewer NA)
#> changed 8,310 values (>99%) of 'Halicryptus' (8310 fewer NA)
#> changed 8,309 values (>99%) of 'Neogobius melanostomus' (8309 fewer NA)
#> changed 8,309 values (>99%) of 'Gobius niger' (8309 fewer NA)
#> changed 8,310 values (>99%) of 'digestive tract' (8310 fewer NA)
#> changed 8,310 values (>99%) of 'Pleuronectidae' (8310 fewer NA)
#> changed 8,310 values (>99%) of 'sprattus sprattus' (8310 fewer NA)
#> changed 8,310 values (>99%) of 'Gasterosteidae' (8310 fewer NA)
#> changed 8,308 values (>99%) of 'litter' (8308 fewer NA)
#> changed 8,289 values (>99%) of 'Mysida' (8289 fewer NA)
#> changed 8,310 values (>99%) of 'Carbon' (8310 fewer NA)
#> changed 8,282 values (>99%) of 'Copepoda' (8282 fewer NA)
#> changed 8,309 values (>99%) of 'Calanoida' (8309 fewer NA)
#> changed 8,310 values (>99%) of 'Belone belone' (8310 fewer NA)
#> changed 8,212 values (99%) of 'Pontoporeiidae' (8212 fewer NA)
#> changed 8,307 values (>99%) of 'Mucus' (8307 fewer NA)
#> changed 8,310 values (>99%) of 'Pectinaria sp.' (8310 fewer NA)
#> changed 8,191 values (99%) of 'Macoma balthica' (8191 fewer NA)
#> changed 8,026 values (97%) of 'remains' (8026 fewer NA)
#> changed 8,258 values (99%) of 'stone' (8258 fewer NA)
#> changed 8,310 values (>99%) of 'carbon' (8310 fewer NA)
#> changed 8,309 values (>99%) of 'Nephtys ciliata' (8309 fewer NA)
#> changed 8,310 values (>99%) of 'Gastrosacus' (8310 fewer NA)
#> changed 8,310 values (>99%) of 'wood' (8310 fewer NA)
#> changed 8,310 values (>99%) of 'Litter/plastics' (8310 fewer NA)
#> changed 8,310 values (>99%) of 'clupeidae' (8310 fewer NA)
#> changed 8,298 values (>99%) of 'Pontoporeidae' (8298 fewer NA)
#> changed 8,304 values (>99%) of 'sand' (8304 fewer NA)
#> changed 8,309 values (>99%) of 'Decapoda' (8309 fewer NA)
#> changed 8,310 values (>99%) of 'Agonus cataphractus' (8310 fewer NA)
#> changed 8,310 values (>99%) of 'clupea harengus' (8310 fewer NA)
#> changed 8,295 values (>99%) of 'Halicryptus spinolusus' (8295 fewer NA)
#> changed 8,156 values (98%) of 'Mytilus edulis' (8156 fewer NA)
#> changed 8,303 values (>99%) of 'priapulida' (8303 fewer NA)
#> changed 8,310 values (>99%) of 'Prapulida' (8310 fewer NA)
#> changed 8,310 values (>99%) of 'Myoxocephalus scorpius' (8310 fewer NA)
head(d_wide)
#> # A tibble: 6 x 113
#> unique_pred_id sprattus_spratt… bylgides_sarsi saduria_entomon clupea_harengus
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 2007_4_101_20… 5.49 0 0 0
#> 2 2007_4_101_20… 0 0.24 0 0
#> 3 2007_4_101_20… 0 0 0.9 0
#> 4 2007_4_101_20… 26.6 0 0 180.
#> 5 2007_4_101_20… 0 0 0 88.5
#> 6 2007_4_101_20… 27.3 0 0 163.
#> # … with 108 more variables: pisces <dbl>, gadus_morhua <dbl>, cumacea <dbl>,
#> # limecola_balthica <dbl>, mysis_mixta <dbl>, corophium_volutator <dbl>,
#> # halicryptus_spinulosus <dbl>, stone <dbl>, pontoporeia_femorata <dbl>,
#> # enchelyopus_cimbrius <dbl>, clupeidae <dbl>, gobiidae <dbl>,
#> # neomysis_integer <dbl>, monoporeia_affinis <dbl>, gammarus_sp <dbl>,
#> # mysidae <dbl>, scoloplos_armiger <dbl>, plastic <dbl>,
#> # priapulus_caudatus <dbl>, bivalvia <dbl>, hyperia_galba <dbl>,
#> # mytilus_sp <dbl>, crangon_crangon <dbl>, idotea_balthica <dbl>,
#> # trachurus_trachurus <dbl>, gasterosteus_aculeatus <dbl>, annelida <dbl>,
#> # algae <dbl>, platichthys_flesus <dbl>, priapulidae <dbl>, waste <dbl>,
#> # praunus_flexuosus <dbl>, unidentified_mass <dbl>,
#> # merlangius_merlangus <dbl>, scales <dbl>, gadidae <dbl>, crustacea <dbl>,
#> # amphipoda <dbl>, spine <dbl>, pleuronectes_platessa <dbl>,
#> # anguilla_anguilla <dbl>, empty <dbl>, mollusca <dbl>,
#> # pungitius_pungitius <dbl>, diastylis_rathkei <dbl>,
#> # terebellides_stroemii <dbl>, palaemon_sp <dbl>, ammodytes_tobianus <dbl>,
#> # na <dbl>, cottidae <dbl>, hediste_diversicolor <dbl>,
#> # palaemon_elegans <dbl>, spinachia_spinachia <dbl>, zoarces_viviparus <dbl>,
#> # ammodytidae <dbl>, caridea <dbl>, amphibalanus_improvisus <dbl>,
#> # myoxocephalus_quadricornis <dbl>, phyllodocida <dbl>, remains <dbl>,
#> # palaemonidae <dbl>, carcinus_maenas <dbl>, gastropoda <dbl>, sand <dbl>,
#> # cerastoderma_glaucum <dbl>, hyperoplus_lanceolatus <dbl>,
#> # mya_arenaria <dbl>, hydrobia_sp <dbl>, wood <dbl>, pleuronectiformes <dbl>,
#> # scophthalmus_maximus <dbl>, polychaeta <dbl>, priapulida <dbl>,
#> # halicryptus <dbl>, neogobius_melanostomus <dbl>, gobius_niger <dbl>,
#> # digestive_tract <dbl>, pleuronectidae <dbl>, sprattus_sprattus_2 <dbl>,
#> # gasterosteidae <dbl>, litter <dbl>, mysida <dbl>, carbon <dbl>,
#> # copepoda <dbl>, calanoida <dbl>, belone_belone <dbl>, pontoporeiidae <dbl>,
#> # mucus <dbl>, pectinaria_sp <dbl>, macoma_balthica <dbl>, remains_2 <dbl>,
#> # stone_2 <dbl>, carbon_2 <dbl>, nephtys_ciliata <dbl>, gastrosacus <dbl>,
#> # wood_2 <dbl>, litter_plastics <dbl>, clupeidae_2 <dbl>,
#> # pontoporeidae <dbl>, sand_2 <dbl>, …
str(d_wide)
#> tibble [8,311 × 113] (S3: tbl_df/tbl/data.frame)
#> $ unique_pred_id : chr [1:8311] "2007_4_101_2007_11_7_101_DAN2_40G5_1_M_490_1028_1_COD" "2007_4_101_2007_11_7_101_DAN2_40G5_118_F_210_82_1_COD" "2007_4_101_2007_11_7_101_DAN2_40G5_45_F_480_784_1_COD" "2007_4_101_2007_11_7_101_DAN2_40G5_65_F_920_7140_1_COD" ...
#> $ sprattus_sprattus : num [1:8311] 5.49 0 0 26.59 0 ...
#> $ bylgides_sarsi : num [1:8311] 0 0.24 0 0 0 0 0 0 0 0 ...
#> $ saduria_entomon : num [1:8311] 0 0 0.9 0 0 0 0 0 0 0 ...
#> $ clupea_harengus : num [1:8311] 0 0 0 179.7 88.5 ...
#> $ pisces : num [1:8311] 0 0 0 20.8 0 ...
#> $ gadus_morhua : num [1:8311] 0 0 0 0 114 ...
#> $ cumacea : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ limecola_balthica : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ mysis_mixta : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ corophium_volutator : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ halicryptus_spinulosus : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ stone : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ pontoporeia_femorata : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ enchelyopus_cimbrius : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ clupeidae : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ gobiidae : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ neomysis_integer : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ monoporeia_affinis : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ gammarus_sp : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ mysidae : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ scoloplos_armiger : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ plastic : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ priapulus_caudatus : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ bivalvia : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ hyperia_galba : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ mytilus_sp : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ crangon_crangon : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ idotea_balthica : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ trachurus_trachurus : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ gasterosteus_aculeatus : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ annelida : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ algae : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ platichthys_flesus : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ priapulidae : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ waste : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ praunus_flexuosus : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ unidentified_mass : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ merlangius_merlangus : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ scales : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ gadidae : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ crustacea : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ amphipoda : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ spine : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ pleuronectes_platessa : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ anguilla_anguilla : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ empty : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ mollusca : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ pungitius_pungitius : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ diastylis_rathkei : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ terebellides_stroemii : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ palaemon_sp : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ ammodytes_tobianus : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ na : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ cottidae : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ hediste_diversicolor : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ palaemon_elegans : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ spinachia_spinachia : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ zoarces_viviparus : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ ammodytidae : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ caridea : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ amphibalanus_improvisus : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ myoxocephalus_quadricornis: num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ phyllodocida : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ remains : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ palaemonidae : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ carcinus_maenas : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ gastropoda : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ sand : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ cerastoderma_glaucum : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ hyperoplus_lanceolatus : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ mya_arenaria : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ hydrobia_sp : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ wood : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ pleuronectiformes : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ scophthalmus_maximus : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ polychaeta : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ priapulida : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ halicryptus : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ neogobius_melanostomus : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ gobius_niger : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ digestive_tract : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ pleuronectidae : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ sprattus_sprattus_2 : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ gasterosteidae : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ litter : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ mysida : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ carbon : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ copepoda : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ calanoida : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ belone_belone : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ pontoporeiidae : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ mucus : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ pectinaria_sp : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ macoma_balthica : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ remains_2 : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ stone_2 : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ carbon_2 : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ nephtys_ciliata : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> [list output truncated]
colnames(d_wide)
#> [1] "unique_pred_id" "sprattus_sprattus"
#> [3] "bylgides_sarsi" "saduria_entomon"
#> [5] "clupea_harengus" "pisces"
#> [7] "gadus_morhua" "cumacea"
#> [9] "limecola_balthica" "mysis_mixta"
#> [11] "corophium_volutator" "halicryptus_spinulosus"
#> [13] "stone" "pontoporeia_femorata"
#> [15] "enchelyopus_cimbrius" "clupeidae"
#> [17] "gobiidae" "neomysis_integer"
#> [19] "monoporeia_affinis" "gammarus_sp"
#> [21] "mysidae" "scoloplos_armiger"
#> [23] "plastic" "priapulus_caudatus"
#> [25] "bivalvia" "hyperia_galba"
#> [27] "mytilus_sp" "crangon_crangon"
#> [29] "idotea_balthica" "trachurus_trachurus"
#> [31] "gasterosteus_aculeatus" "annelida"
#> [33] "algae" "platichthys_flesus"
#> [35] "priapulidae" "waste"
#> [37] "praunus_flexuosus" "unidentified_mass"
#> [39] "merlangius_merlangus" "scales"
#> [41] "gadidae" "crustacea"
#> [43] "amphipoda" "spine"
#> [45] "pleuronectes_platessa" "anguilla_anguilla"
#> [47] "empty" "mollusca"
#> [49] "pungitius_pungitius" "diastylis_rathkei"
#> [51] "terebellides_stroemii" "palaemon_sp"
#> [53] "ammodytes_tobianus" "na"
#> [55] "cottidae" "hediste_diversicolor"
#> [57] "palaemon_elegans" "spinachia_spinachia"
#> [59] "zoarces_viviparus" "ammodytidae"
#> [61] "caridea" "amphibalanus_improvisus"
#> [63] "myoxocephalus_quadricornis" "phyllodocida"
#> [65] "remains" "palaemonidae"
#> [67] "carcinus_maenas" "gastropoda"
#> [69] "sand" "cerastoderma_glaucum"
#> [71] "hyperoplus_lanceolatus" "mya_arenaria"
#> [73] "hydrobia_sp" "wood"
#> [75] "pleuronectiformes" "scophthalmus_maximus"
#> [77] "polychaeta" "priapulida"
#> [79] "halicryptus" "neogobius_melanostomus"
#> [81] "gobius_niger" "digestive_tract"
#> [83] "pleuronectidae" "sprattus_sprattus_2"
#> [85] "gasterosteidae" "litter"
#> [87] "mysida" "carbon"
#> [89] "copepoda" "calanoida"
#> [91] "belone_belone" "pontoporeiidae"
#> [93] "mucus" "pectinaria_sp"
#> [95] "macoma_balthica" "remains_2"
#> [97] "stone_2" "carbon_2"
#> [99] "nephtys_ciliata" "gastrosacus"
#> [101] "wood_2" "litter_plastics"
#> [103] "clupeidae_2" "pontoporeidae"
#> [105] "sand_2" "decapoda"
#> [107] "agonus_cataphractus" "clupea_harengus_2"
#> [109] "halicryptus_spinolusus" "mytilus_edulis"
#> [111] "priapulida_2" "prapulida"
#> [113] "myoxocephalus_scorpius"
# Now make some calculations and aggregate some taxonomic level
d_wide2 <- d_wide %>%
mutate(amphipoda_tot = bylgides_sarsi + hyperia_galba + gammarus_sp + monoporeia_affinis +
corophium_volutator + amphipoda,
bivalvia_tot = bivalvia + mytilus_sp + cerastoderma_glaucum + mya_arenaria + macoma_balthica +
mytilus_edulis,
clupeidae_tot = clupeidae + clupea_harengus_2 + sprattus_sprattus_2 + clupeidae_2,
sprattus_sprattus_tot = sprattus_sprattus + sprattus_sprattus_2,
clupea_harengus_tot = clupea_harengus + clupea_harengus_2,
#gadus_morhua_tot = gadus_morhua,
gadiformes_tot = gadidae + merlangius_merlangus,
gobiidae_tot = gobiidae,
mysidae_tot = mysidae + neomysis_integer + mysis_mixta + mysida + gastrosacus,
non_bio_tot = stone + plastic + sand + wood + litter + carbon + stone_2 + carbon_2 + wood_2 +
litter_plastics + sand_2,
other_crustacea_tot = pontoporeia_femorata + crangon_crangon + idotea_balthica +
praunus_flexuosus + crustacea + diastylis_rathkei + palaemon_sp + palaemon_elegans + caridea +
amphibalanus_improvisus + palaemonidae + carcinus_maenas + copepoda + calanoida + pontoporeiidae + decapoda,
other_tot = halicryptus_spinulosus + priapulus_caudatus + annelida + algae + priapulidae + waste +
unidentified_mass + spine + empty + mollusca + na + remains + gastropoda + hydrobia_sp +
priapulida + halicryptus + digestive_tract + mucus + remains_2 + pontoporeidae +
halicryptus_spinolusus + priapulida_2 + prapulida,
other_pisces_tot = enchelyopus_cimbrius + trachurus_trachurus + gasterosteus_aculeatus +
scales + pleuronectes_platessa + anguilla_anguilla + pungitius_pungitius + ammodytes_tobianus +
cottidae + spinachia_spinachia + zoarces_viviparus + ammodytidae + myoxocephalus_quadricornis +
hyperoplus_lanceolatus + pleuronectiformes + scophthalmus_maximus + neogobius_melanostomus + gobius_niger +
pleuronectidae + gasterosteidae + belone_belone + agonus_cataphractus + myoxocephalus_scorpius,
platichthys_flesus_tot = platichthys_flesus,
polychaeta_tot = bylgides_sarsi + scoloplos_armiger + terebellides_stroemii + hediste_diversicolor +
phyllodocida + polychaeta + pectinaria_sp + nephtys_ciliata,
saduria_entomon_tot = saduria_entomon
)
#> mutate: new variable 'amphipoda_tot' (double) with 397 unique values and 0% NA
#> new variable 'bivalvia_tot' (double) with 257 unique values and 0% NA
#> new variable 'clupeidae_tot' (double) with 285 unique values and 0% NA
#> new variable 'sprattus_sprattus_tot' (double) with 986 unique values and 0% NA
#> new variable 'clupea_harengus_tot' (double) with 646 unique values and 0% NA
#> new variable 'gadiformes_tot' (double) with 7 unique values and 0% NA
#> new variable 'gobiidae_tot' (double) with 163 unique values and 0% NA
#> new variable 'mysidae_tot' (double) with 358 unique values and 0% NA
#> new variable 'non_bio_tot' (double) with 61 unique values and 0% NA
#> new variable 'other_crustacea_tot' (double) with 285 unique values and 0% NA
#> new variable 'other_tot' (double) with 299 unique values and 0% NA
#> new variable 'other_pisces_tot' (double) with 253 unique values and 0% NA
#> new variable 'platichthys_flesus_tot' (double) with 16 unique values and 0% NA
#> new variable 'polychaeta_tot' (double) with 335 unique values and 0% NA
#> new variable 'saduria_entomon_tot' (double) with 626 unique values and 0% NA
length(unique(d$Unique.pred.id))
#> [1] 9432
str(d_wide2) # Check why they do not match in length / nrows
#> tibble [8,311 × 128] (S3: tbl_df/tbl/data.frame)
#> $ unique_pred_id : chr [1:8311] "2007_4_101_2007_11_7_101_DAN2_40G5_1_M_490_1028_1_COD" "2007_4_101_2007_11_7_101_DAN2_40G5_118_F_210_82_1_COD" "2007_4_101_2007_11_7_101_DAN2_40G5_45_F_480_784_1_COD" "2007_4_101_2007_11_7_101_DAN2_40G5_65_F_920_7140_1_COD" ...
#> $ sprattus_sprattus : num [1:8311] 5.49 0 0 26.59 0 ...
#> $ bylgides_sarsi : num [1:8311] 0 0.24 0 0 0 0 0 0 0 0 ...
#> $ saduria_entomon : num [1:8311] 0 0 0.9 0 0 0 0 0 0 0 ...
#> $ clupea_harengus : num [1:8311] 0 0 0 179.7 88.5 ...
#> $ pisces : num [1:8311] 0 0 0 20.8 0 ...
#> $ gadus_morhua : num [1:8311] 0 0 0 0 114 ...
#> $ cumacea : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ limecola_balthica : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ mysis_mixta : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ corophium_volutator : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ halicryptus_spinulosus : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ stone : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ pontoporeia_femorata : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ enchelyopus_cimbrius : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ clupeidae : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ gobiidae : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ neomysis_integer : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ monoporeia_affinis : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ gammarus_sp : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ mysidae : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ scoloplos_armiger : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ plastic : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ priapulus_caudatus : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ bivalvia : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ hyperia_galba : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ mytilus_sp : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ crangon_crangon : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ idotea_balthica : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ trachurus_trachurus : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ gasterosteus_aculeatus : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ annelida : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ algae : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ platichthys_flesus : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ priapulidae : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ waste : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ praunus_flexuosus : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ unidentified_mass : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ merlangius_merlangus : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ scales : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ gadidae : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ crustacea : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ amphipoda : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ spine : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ pleuronectes_platessa : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ anguilla_anguilla : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ empty : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ mollusca : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ pungitius_pungitius : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ diastylis_rathkei : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ terebellides_stroemii : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ palaemon_sp : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ ammodytes_tobianus : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ na : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ cottidae : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ hediste_diversicolor : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ palaemon_elegans : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ spinachia_spinachia : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ zoarces_viviparus : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ ammodytidae : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ caridea : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ amphibalanus_improvisus : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ myoxocephalus_quadricornis: num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ phyllodocida : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ remains : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ palaemonidae : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ carcinus_maenas : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ gastropoda : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ sand : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ cerastoderma_glaucum : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ hyperoplus_lanceolatus : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ mya_arenaria : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ hydrobia_sp : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ wood : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ pleuronectiformes : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ scophthalmus_maximus : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ polychaeta : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ priapulida : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ halicryptus : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ neogobius_melanostomus : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ gobius_niger : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ digestive_tract : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ pleuronectidae : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ sprattus_sprattus_2 : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ gasterosteidae : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ litter : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ mysida : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ carbon : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ copepoda : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ calanoida : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ belone_belone : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ pontoporeiidae : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ mucus : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ pectinaria_sp : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ macoma_balthica : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ remains_2 : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ stone_2 : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ carbon_2 : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> $ nephtys_ciliata : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#> [list output truncated]
# Select only columns aggregated columns
d_wide3 <- d_wide2 %>% dplyr::select(c(1, 114:128))
colnames(d_wide3)
#> [1] "unique_pred_id" "amphipoda_tot" "bivalvia_tot"
#> [4] "clupeidae_tot" "sprattus_sprattus_tot" "clupea_harengus_tot"
#> [7] "gadiformes_tot" "gobiidae_tot" "mysidae_tot"
#> [10] "non_bio_tot" "other_crustacea_tot" "other_tot"
#> [13] "other_pisces_tot" "platichthys_flesus_tot" "polychaeta_tot"
#> [16] "saduria_entomon_tot"
# Add back in other information
d_sub <- d %>%
dplyr::select(Year, Quarter, Cruise, Predator.code, X, Y, Lat, Long, Ices.rect, Pred.size.mm,
Pred.weight.g, Unique.pred.id, source) %>%
rename("Predator.spec" = "Predator.code") %>%
distinct(Unique.pred.id, .keep_all = TRUE) %>%
janitor::clean_names()
#> rename: renamed one variable (Predator.spec)
#> distinct: removed 12,299 rows (57%), 9,432 rows remaining
d_wide4 <- left_join(d_wide3, d_sub) # Why missing 1,121 IDs?
#> Joining, by = "unique_pred_id"
#> left_join: added 12 columns (year, quarter, cruise, predator_spec, x, …)
#> > rows only in x 0
#> > rows only in y (1,121)
#> > matched rows 8,311
#> > =======
#> > rows total 8,311
#d_wide4 %>% group_by(unique_pred_id) %>% mutate(n = n()) %>% ungroup() %>% distinct(n)
colnames(d_wide4)
#> [1] "unique_pred_id" "amphipoda_tot" "bivalvia_tot"
#> [4] "clupeidae_tot" "sprattus_sprattus_tot" "clupea_harengus_tot"
#> [7] "gadiformes_tot" "gobiidae_tot" "mysidae_tot"
#> [10] "non_bio_tot" "other_crustacea_tot" "other_tot"
#> [13] "other_pisces_tot" "platichthys_flesus_tot" "polychaeta_tot"
#> [16] "saduria_entomon_tot" "year" "quarter"
#> [19] "cruise" "predator_spec" "x"
#> [22] "y" "lat" "long"
#> [25] "ices_rect" "pred_size_mm" "pred_weight_g"
#> [28] "source"
head(data.frame(d_wide4))
#> unique_pred_id amphipoda_tot
#> 1 2007_4_101_2007_11_7_101_DAN2_40G5_1_M_490_1028_1_COD 0.00
#> 2 2007_4_101_2007_11_7_101_DAN2_40G5_118_F_210_82_1_COD 0.24
#> 3 2007_4_101_2007_11_7_101_DAN2_40G5_45_F_480_784_1_COD 0.00
#> 4 2007_4_101_2007_11_7_101_DAN2_40G5_65_F_920_7140_1_COD 0.00
#> 5 2007_4_101_2007_11_7_101_DAN2_40G5_66_M_740_3812_1_COD 0.00
#> 6 2007_4_101_2007_11_7_101_DAN2_40G5_70_F_590_2198_1_COD 0.00
#> bivalvia_tot clupeidae_tot sprattus_sprattus_tot clupea_harengus_tot
#> 1 0 0 5.49 0.00
#> 2 0 0 0.00 0.00
#> 3 0 0 0.00 0.00
#> 4 0 0 26.59 179.66
#> 5 0 0 0.00 88.51
#> 6 0 0 27.32 163.47
#> gadiformes_tot gobiidae_tot mysidae_tot non_bio_tot other_crustacea_tot
#> 1 0 0 0 0 0
#> 2 0 0 0 0 0
#> 3 0 0 0 0 0
#> 4 0 0 0 0 0
#> 5 0 0 0 0 0
#> 6 0 0 0 0 0
#> other_tot other_pisces_tot platichthys_flesus_tot polychaeta_tot
#> 1 0 0 0 0.00
#> 2 0 0 0 0.24
#> 3 0 0 0 0.00
#> 4 0 0 0 0.00
#> 5 0 0 0 0.00
#> 6 0 0 0 0.00
#> saduria_entomon_tot year quarter cruise predator_spec x y
#> 1 0.0 2007 4 BITS-2 COD 520.8793 6176.115
#> 2 0.0 2007 4 BITS-2 COD 520.8793 6176.115
#> 3 0.9 2007 4 BITS-2 COD 520.8793 6176.115
#> 4 0.0 2007 4 BITS-2 COD 520.8793 6176.115
#> 5 0.0 2007 4 BITS-2 COD 520.8793 6176.115
#> 6 0.0 2007 4 BITS-2 COD 520.8793 6176.115
#> lat long ices_rect pred_size_mm pred_weight_g source
#> 1 55.73032 15.33247 40G5 490 1028 Kevin's MSC
#> 2 55.73032 15.33247 40G5 210 82 Kevin's MSC
#> 3 55.73032 15.33247 40G5 480 784 Kevin's MSC
#> 4 55.73032 15.33247 40G5 920 7140 Kevin's MSC
#> 5 55.73032 15.33247 40G5 740 3812 Kevin's MSC
#> 6 55.73032 15.33247 40G5 590 2198 Kevin's MSC
colnames(d_wide4)
#> [1] "unique_pred_id" "amphipoda_tot" "bivalvia_tot"
#> [4] "clupeidae_tot" "sprattus_sprattus_tot" "clupea_harengus_tot"
#> [7] "gadiformes_tot" "gobiidae_tot" "mysidae_tot"
#> [10] "non_bio_tot" "other_crustacea_tot" "other_tot"
#> [13] "other_pisces_tot" "platichthys_flesus_tot" "polychaeta_tot"
#> [16] "saduria_entomon_tot" "year" "quarter"
#> [19] "cruise" "predator_spec" "x"
#> [22] "y" "lat" "long"
#> [25] "ices_rect" "pred_size_mm" "pred_weight_g"
#> [28] "source"
d_wide4[1, 2]
#> # A tibble: 1 x 1
#> amphipoda_tot
#> <dbl>
#> 1 0
d_wide4[1, 16]
#> # A tibble: 1 x 1
#> saduria_entomon_tot
#> <dbl>
#> 1 0
d_wide4 <- d_wide4 %>% mutate(tot_prey_biom = rowSums(.[2:16]))
#> mutate: new variable 'tot_prey_biom' (double) with 2,467 unique values and 0% NA
# Now, group by survey and Quarter, select only the main groups (as Haase) and see if it resemlbes the
# previous results.
fle <- d_wide4 %>%
filter(predator_spec == "FLE") %>%
mutate(pred_cm = pred_size_mm/10,
pred_cm_class = cut(pred_cm, breaks = c(0, 9, 20, 30, 200), right = TRUE))
#> filter: removed 6,054 rows (73%), 2,257 rows remaining
#> mutate: new variable 'pred_cm' (double) with 113 unique values and <1% NA
#> new variable 'pred_cm_class' (factor) with 5 unique values and <1% NA
head(data.frame(fle), 20)
#> unique_pred_id amphipoda_tot bivalvia_tot clupeidae_tot
#> 1 2015_2_1010_4_FLE 0 2.60 0
#> 2 2015_2_1010_402_FLE 0 0.00 0
#> 3 2015_2_108_3_FLE 0 0.00 0
#> 4 2015_2_202_403_FLE 0 5.46 0
#> 5 2015_2_203_1_FLE 0 1.74 0
#> 6 2015_2_204_2_FLE 0 0.04 0
#> 7 2015_2_501_10_FLE 0 0.99 0
#> 8 2015_2_501_60_FLE 0 0.00 0
#> 9 2015_2_501_8_FLE 0 0.00 0
#> 10 2015_2_501_9_FLE 0 0.00 0
#> 11 2015_2_504_5_FLE 0 0.00 0
#> 12 2015_2_504_58_FLE 0 0.00 0
#> 13 2015_2_504_6_FLE 0 0.00 0
#> 14 2015_2_504_7_FLE 0 0.00 0
#> 15 2015_2_506_11_FLE 0 0.00 0
#> 16 2015_2_506_12_FLE 0 0.00 0
#> 17 2015_2_506_13_FLE 0 0.00 0
#> 18 2015_2_506_14_FLE 0 0.00 0
#> 19 2015_2_506_15_FLE 0 0.00 0
#> 20 2015_2_506_16_FLE 0 0.00 0
#> sprattus_sprattus_tot clupea_harengus_tot gadiformes_tot gobiidae_tot
#> 1 0 0 0 0
#> 2 0 0 0 0
#> 3 0 0 0 0
#> 4 0 0 0 0
#> 5 0 0 0 0
#> 6 0 0 0 0
#> 7 0 0 0 0
#> 8 0 0 0 0
#> 9 0 0 0 0
#> 10 0 0 0 0
#> 11 0 0 0 0
#> 12 0 0 0 0
#> 13 0 0 0 0
#> 14 0 0 0 0
#> 15 0 0 0 0
#> 16 0 0 0 0
#> 17 0 0 0 0
#> 18 0 0 0 0
#> 19 0 0 0 0
#> 20 0 0 0 0
#> mysidae_tot non_bio_tot other_crustacea_tot other_tot other_pisces_tot
#> 1 0 0.19 0.00 0.00 2.93
#> 2 0 0.00 0.00 0.00 0.00
#> 3 0 0.00 0.00 0.00 0.00
#> 4 0 0.13 0.00 0.00 0.00
#> 5 0 0.00 0.00 0.00 0.00
#> 6 0 0.00 0.00 0.00 0.00
#> 7 0 0.00 0.00 0.00 0.00
#> 8 0 0.00 0.00 0.00 0.00
#> 9 0 0.00 0.01 0.03 0.00
#> 10 0 0.00 0.00 0.00 0.00
#> 11 0 0.00 0.02 0.00 0.00
#> 12 0 0.00 0.00 0.02 0.00
#> 13 0 0.00 0.00 0.00 0.00
#> 14 0 0.00 0.00 0.00 0.00
#> 15 0 0.00 0.00 0.00 0.00
#> 16 0 0.00 0.10 0.00 0.00
#> 17 0 0.00 0.00 0.00 0.00
#> 18 0 0.00 0.00 0.00 0.00
#> 19 0 0.00 0.00 0.00 0.00
#> 20 0 0.00 0.00 0.00 0.00
#> platichthys_flesus_tot polychaeta_tot saduria_entomon_tot year quarter
#> 1 0 0 0.00 2015 2
#> 2 0 0 0.00 2015 2
#> 3 0 0 5.61 2015 2
#> 4 0 0 1.02 2015 2
#> 5 0 0 3.10 2015 2
#> 6 0 0 2.98 2015 2
#> 7 0 0 0.00 2015 2
#> 8 0 0 0.00 2015 2
#> 9 0 0 0.00 2015 2
#> 10 0 0 0.00 2015 2
#> 11 0 0 0.08 2015 2
#> 12 0 0 0.00 2015 2
#> 13 0 0 0.00 2015 2
#> 14 0 0 0.00 2015 2
#> 15 0 0 0.00 2015 2
#> 16 0 0 0.00 2015 2
#> 17 0 0 0.00 2015 2
#> 18 0 0 0.00 2015 2
#> 19 0 0 0.00 2015 2
#> 20 0 0 0.00 2015 2
#> cruise predator_spec x y lat long ices_rect pred_size_mm
#> 1 INSPIRE FLE 486.2897 6209.440 56.03 14.78 41G4 270
#> 2 INSPIRE FLE 486.2897 6209.440 56.03 14.78 41G4 304
#> 3 INSPIRE FLE 483.1607 6206.112 56.00 14.73 41G4 350
#> 4 INSPIRE FLE 485.6369 6200.539 55.95 14.77 40G4 302
#> 5 INSPIRE FLE 486.9129 6209.438 56.03 14.79 41G4 299
#> 6 INSPIRE FLE 485.6517 6204.990 55.99 14.77 40G4 346
#> 7 INSPIRE FLE 478.0250 6177.198 55.74 14.65 40G4 314
#> 8 INSPIRE FLE 478.0250 6177.198 55.74 14.65 40G4 330
#> 9 INSPIRE FLE 478.0250 6177.198 55.74 14.65 40G4 202
#> 10 INSPIRE FLE 478.0250 6177.198 55.74 14.65 40G4 298
#> 11 INSPIRE FLE 477.9969 6171.634 55.69 14.65 40G4 352
#> 12 INSPIRE FLE 477.9969 6171.634 55.69 14.65 40G4 313
#> 13 INSPIRE FLE 477.9969 6171.634 55.69 14.65 40G4 328
#> 14 INSPIRE FLE 477.9969 6171.634 55.69 14.65 40G4 250
#> 15 INSPIRE FLE 479.2754 6176.079 55.73 14.67 40G4 299
#> 16 INSPIRE FLE 479.2754 6176.079 55.73 14.67 40G4 280
#> 17 INSPIRE FLE 479.2754 6176.079 55.73 14.67 40G4 349
#> 18 INSPIRE FLE 479.2754 6176.079 55.73 14.67 40G4 335
#> 19 INSPIRE FLE 479.2754 6176.079 55.73 14.67 40G4 340
#> 20 INSPIRE FLE 479.2754 6176.079 55.73 14.67 40G4 337
#> pred_weight_g source tot_prey_biom pred_cm pred_cm_class
#> 1 236 Kevin's MSC 5.72 27.0 (20,30]
#> 2 476 Kevin's MSC 0.00 30.4 (30,200]
#> 3 334 Kevin's MSC 5.61 35.0 (30,200]
#> 4 264 Kevin's MSC 6.61 30.2 (30,200]
#> 5 246 Kevin's MSC 4.84 29.9 (20,30]
#> 6 346 Kevin's MSC 3.02 34.6 (30,200]
#> 7 305 Kevin's MSC 0.99 31.4 (30,200]
#> 8 548 Kevin's MSC 0.00 33.0 (30,200]
#> 9 89 Kevin's MSC 0.04 20.2 (20,30]
#> 10 251 Kevin's MSC 0.00 29.8 (20,30]
#> 11 325 Kevin's MSC 0.10 35.2 (30,200]
#> 12 319 Kevin's MSC 0.02 31.3 (30,200]
#> 13 367 Kevin's MSC 0.00 32.8 (30,200]
#> 14 172 Kevin's MSC 0.00 25.0 (20,30]
#> 15 196 Kevin's MSC 0.00 29.9 (20,30]
#> 16 184 Kevin's MSC 0.10 28.0 (20,30]
#> 17 389 Kevin's MSC 0.00 34.9 (30,200]
#> 18 368 Kevin's MSC 0.00 33.5 (30,200]
#> 19 371 Kevin's MSC 0.00 34.0 (30,200]
#> 20 380 Kevin's MSC 0.00 33.7 (30,200]
unique(fle$pred_cm_class)
#> [1] (20,30] (30,200] (9,20] (0,9] <NA>
#> Levels: (0,9] (9,20] (20,30] (30,200]
cod <- d_wide4 %>%
filter(predator_spec == "COD") %>%
mutate(pred_cm = pred_size_mm/10,
pred_cm_class = cut(pred_cm, breaks = c(0, 6, 20, 30, 40, 50, 200), right = TRUE))
#> filter: removed 2,257 rows (27%), 6,054 rows remaining
#> mutate: new variable 'pred_cm' (double) with 222 unique values and <1% NA
#> new variable 'pred_cm_class' (factor) with 7 unique values and <1% NA
head(data.frame(cod), 20)
#> unique_pred_id amphipoda_tot
#> 1 2007_4_101_2007_11_7_101_DAN2_40G5_1_M_490_1028_1_COD 0.00
#> 2 2007_4_101_2007_11_7_101_DAN2_40G5_118_F_210_82_1_COD 0.24
#> 3 2007_4_101_2007_11_7_101_DAN2_40G5_45_F_480_784_1_COD 0.00
#> 4 2007_4_101_2007_11_7_101_DAN2_40G5_65_F_920_7140_1_COD 0.00
#> 5 2007_4_101_2007_11_7_101_DAN2_40G5_66_M_740_3812_1_COD 0.00
#> 6 2007_4_101_2007_11_7_101_DAN2_40G5_70_F_590_2198_1_COD 0.00
#> 7 2007_4_101_2007_11_7_101_DAN2_40G5_72_F_540_1312_1_COD 0.00
#> 8 2007_4_101_2007_11_7_101_DAN2_40G5_73_F_560_1618_1_COD 0.00
#> 9 2007_4_101_2007_11_7_101_DAN2_40G5_75_F_780_5026_1_COD 0.00
#> 10 2007_4_101_2007_11_7_101_DAN2_40G5_76_M_590_1868_1_COD 0.00
#> 11 2007_4_101_2007_11_7_101_DAN2_40G5_77_F_590_1320_1_COD 0.00
#> 12 2007_4_101_2007_11_7_101_DAN2_40G5_81_M_630_2002_1_COD 0.00
#> 13 2007_4_101_2007_11_7_101_DAN2_40G5_82_F_680_2934_1_COD 0.00
#> 14 2007_4_101_2007_11_7_101_DAN2_40G5_84_F_580_1382_1_COD 0.00
#> 15 2007_4_101_2007_11_7_101_DAN2_40G5_86_F_580_1360_1_COD 0.00
#> 16 2007_4_101_2007_11_7_101_DAN2_40G5_87_F_540_1380_1_COD 0.00
#> 17 2007_4_101_2007_11_7_101_DAN2_40G5_94_M_540_1536_1_COD 0.00
#> 18 2007_4_101_2007_11_7_101_DAN2_40G5_95_M_590_1714_1_COD 0.00
#> 19 2007_4_103_2007_11_7_103_DAN2_40G5_1_M_740_3350_1_COD 0.00
#> 20 2007_4_103_2007_11_7_103_DAN2_40G5_122_F_650_2588_1_COD 0.00
#> bivalvia_tot clupeidae_tot sprattus_sprattus_tot clupea_harengus_tot
#> 1 0 0 5.49 0.00
#> 2 0 0 0.00 0.00
#> 3 0 0 0.00 0.00
#> 4 0 0 26.59 179.66
#> 5 0 0 0.00 88.51
#> 6 0 0 27.32 163.47
#> 7 0 0 4.93 0.00
#> 8 0 0 14.67 0.00
#> 9 0 0 51.84 264.95
#> 10 0 0 1.72 0.00
#> 11 0 0 18.24 0.00
#> 12 0 0 12.32 51.80
#> 13 0 0 31.46 169.98
#> 14 0 0 0.00 0.00
#> 15 0 0 19.95 0.00
#> 16 0 0 0.00 0.00
#> 17 0 0 12.68 59.26
#> 18 0 0 24.77 75.62
#> 19 0 0 0.00 0.00
#> 20 0 0 0.00 0.00
#> gadiformes_tot gobiidae_tot mysidae_tot non_bio_tot other_crustacea_tot
#> 1 0 0 0 0 0
#> 2 0 0 0 0 0
#> 3 0 0 0 0 0
#> 4 0 0 0 0 0
#> 5 0 0 0 0 0
#> 6 0 0 0 0 0
#> 7 0 0 0 0 0
#> 8 0 0 0 0 0
#> 9 0 0 0 0 0
#> 10 0 0 0 0 0
#> 11 0 0 0 0 0
#> 12 0 0 0 0 0
#> 13 0 0 0 0 0
#> 14 0 0 0 0 0
#> 15 0 0 0 0 0
#> 16 0 0 0 0 0
#> 17 0 0 0 0 0
#> 18 0 0 0 0 0
#> 19 0 0 0 0 0
#> 20 0 0 0 0 0
#> other_tot other_pisces_tot platichthys_flesus_tot polychaeta_tot
#> 1 0 0 0 0.00
#> 2 0 0 0 0.24
#> 3 0 0 0 0.00
#> 4 0 0 0 0.00
#> 5 0 0 0 0.00
#> 6 0 0 0 0.00
#> 7 0 0 0 0.00
#> 8 0 0 0 0.00
#> 9 0 0 0 0.00
#> 10 0 0 0 0.00
#> 11 0 0 0 0.00
#> 12 0 0 0 0.00
#> 13 0 0 0 0.00
#> 14 0 0 0 0.00
#> 15 0 0 0 0.00
#> 16 0 0 0 0.00
#> 17 0 0 0 0.00
#> 18 0 0 0 0.00
#> 19 0 0 0 0.00
#> 20 0 0 0 0.00
#> saduria_entomon_tot year quarter cruise predator_spec x y
#> 1 0.000 2007 4 BITS-2 COD 520.8793 6176.115
#> 2 0.000 2007 4 BITS-2 COD 520.8793 6176.115
#> 3 0.900 2007 4 BITS-2 COD 520.8793 6176.115
#> 4 0.000 2007 4 BITS-2 COD 520.8793 6176.115
#> 5 0.000 2007 4 BITS-2 COD 520.8793 6176.115
#> 6 0.000 2007 4 BITS-2 COD 520.8793 6176.115
#> 7 0.000 2007 4 BITS-2 COD 520.8793 6176.115
#> 8 0.000 2007 4 BITS-2 COD 520.8793 6176.115
#> 9 0.000 2007 4 BITS-2 COD 520.8793 6176.115
#> 10 0.000 2007 4 BITS-2 COD 520.8793 6176.115
#> 11 0.950 2007 4 BITS-2 COD 520.8793 6176.115
#> 12 0.000 2007 4 BITS-2 COD 520.8793 6176.115
#> 13 0.000 2007 4 BITS-2 COD 520.8793 6176.115
#> 14 0.000 2007 4 BITS-2 COD 520.8793 6176.115
#> 15 0.000 2007 4 BITS-2 COD 520.8793 6176.115
#> 16 14.590 2007 4 BITS-2 COD 520.8793 6176.115
#> 17 8.691 2007 4 BITS-2 COD 520.8793 6176.115
#> 18 0.000 2007 4 BITS-2 COD 520.8793 6176.115
#> 19 58.260 2007 4 BITS-2 COD 524.3767 6182.403
#> 20 25.120 2007 4 BITS-2 COD 524.3767 6182.403
#> lat long ices_rect pred_size_mm pred_weight_g source
#> 1 55.73032 15.33247 40G5 490 1028 Kevin's MSC
#> 2 55.73032 15.33247 40G5 210 82 Kevin's MSC
#> 3 55.73032 15.33247 40G5 480 784 Kevin's MSC
#> 4 55.73032 15.33247 40G5 920 7140 Kevin's MSC
#> 5 55.73032 15.33247 40G5 740 3812 Kevin's MSC
#> 6 55.73032 15.33247 40G5 590 2198 Kevin's MSC
#> 7 55.73032 15.33247 40G5 540 1312 Kevin's MSC
#> 8 55.73032 15.33247 40G5 560 1618 Kevin's MSC
#> 9 55.73032 15.33247 40G5 780 5026 Kevin's MSC
#> 10 55.73032 15.33247 40G5 590 1868 Kevin's MSC
#> 11 55.73032 15.33247 40G5 590 1320 Kevin's MSC
#> 12 55.73032 15.33247 40G5 630 2002 Kevin's MSC
#> 13 55.73032 15.33247 40G5 680 2934 Kevin's MSC
#> 14 55.73032 15.33247 40G5 580 1382 Kevin's MSC
#> 15 55.73032 15.33247 40G5 580 1360 Kevin's MSC
#> 16 55.73032 15.33247 40G5 540 1380 Kevin's MSC
#> 17 55.73032 15.33247 40G5 540 1536 Kevin's MSC
#> 18 55.73032 15.33247 40G5 590 1714 Kevin's MSC
#> 19 55.78665 15.38872 40G5 740 3350 Kevin's MSC
#> 20 55.78665 15.38872 40G5 650 2588 Kevin's MSC
#> tot_prey_biom pred_cm pred_cm_class
#> 1 5.490 49 (40,50]
#> 2 0.480 21 (20,30]
#> 3 0.900 48 (40,50]
#> 4 206.250 92 (50,200]
#> 5 88.510 74 (50,200]
#> 6 190.790 59 (50,200]
#> 7 4.930 54 (50,200]
#> 8 14.670 56 (50,200]
#> 9 316.790 78 (50,200]
#> 10 1.720 59 (50,200]
#> 11 19.190 59 (50,200]
#> 12 64.120 63 (50,200]
#> 13 201.440 68 (50,200]
#> 14 0.000 58 (50,200]
#> 15 19.950 58 (50,200]
#> 16 14.590 54 (50,200]
#> 17 80.631 54 (50,200]
#> 18 100.390 59 (50,200]
#> 19 58.260 74 (50,200]
#> 20 25.120 65 (50,200]
unique(cod$pred_cm_class)
#> [1] (40,50] (20,30] (50,200] (6,20] (30,40] (0,6] <NA>
#> Levels: (0,6] (6,20] (20,30] (30,40] (40,50] (50,200]Reproduce Haase figures with biomass proportions
# Try to reproduce Haase et al!
# Aggregate across predator individuals and convert to long data for plotting
cod2 <- cod %>%
drop_na(pred_cm_class) %>% # This is because not all predators have length
filter(cruise %in% c("BITS", "BITS-1", "BITS-2")) %>%
filter(year %in% c(2015, 2016, 2017)) %>%
filter(source == "Kevin's MSC") %>%
filter(!pred_cm_class == "(0,6]") %>%
group_by(quarter, pred_cm_class) %>%
summarise(sprattus_sprattus_tot_all = sum(sprattus_sprattus_tot),
saduria_entomon_tot_all = sum(saduria_entomon_tot),
mysidae_tot_all = sum(mysidae_tot),
clupea_harengus_tot_all = sum(clupea_harengus_tot),
sprattus_sprattus_tot_all = sum(sprattus_sprattus_tot),
gobiidae_tot_all = sum(gobiidae_tot),
gadiformes_tot_all = sum(gadiformes_tot),
other_pisces_tot_all = sum(other_pisces_tot),
polychaeta_tot_all = sum(polychaeta_tot)) %>%
ungroup() %>%
pivot_longer(3:10) %>%
group_by(quarter, pred_cm_class) %>%
mutate(percent = value/sum(value)) %>%
ungroup()
#> drop_na: removed 24 rows (<1%), 6,030 rows remaining
#> filter: removed 749 rows (12%), 5,281 rows remaining
#> filter: removed 3,158 rows (60%), 2,123 rows remaining
#> filter: removed 283 rows (13%), 1,840 rows remaining
#> filter: removed 16 rows (1%), 1,824 rows remaining
#> group_by: 2 grouping variables (quarter, pred_cm_class)
#> summarise: now 10 rows and 10 columns, one group variable remaining (quarter)
#> ungroup: no grouping variables
#> pivot_longer: reorganized (sprattus_sprattus_tot_all, saduria_entomon_tot_all, mysidae_tot_all, clupea_harengus_tot_all, gobiidae_tot_all, …) into (name, value) [was 10x10, now 80x4]
#> group_by: 2 grouping variables (quarter, pred_cm_class)
#> mutate (grouped): new variable 'percent' (double) with 65 unique values and 0% NA
#> ungroup: no grouping variables
# So much herring for small cod?
# A 20 cm cod weighs approx 100 g. Can it have 20g of herring in the stomach?
#small_cod <-
d %>%
filter(Pred.size.mm < 200 & Prey.sp.code %in% c("clupea harengus", "Clupea harengus")) %>%
dplyr::select(Prey.weight, Prey.sp.code, Unique.pred.id, Pred.size.mm)
#> filter: removed 21,728 rows (>99%), 3 rows remaining
#> # A tibble: 3 x 4
#> Prey.weight Prey.sp.code Unique.pred.id Pred.size.mm
#> <dbl> <chr> <chr> <dbl>
#> 1 20.0 Clupea hareng… 2016_1_73_735_COD 180
#> 2 4.48 Clupea hareng… 2007_4_126_2007_11_8_126_DAN2_39G4_1_… 160
#> 3 9.18 Clupea hareng… 2007_4_126_2007_11_8_126_DAN2_39G4_4_… 170
d_wide4 %>% filter(unique_pred_id == "2016_1_73_735_COD") %>% as.data.frame()
#> filter: removed 8,310 rows (>99%), one row remaining
#> unique_pred_id amphipoda_tot bivalvia_tot clupeidae_tot
#> 1 2016_1_73_735_COD 0.05 0 0
#> sprattus_sprattus_tot clupea_harengus_tot gadiformes_tot gobiidae_tot
#> 1 0 19.99 0 0
#> mysidae_tot non_bio_tot other_crustacea_tot other_tot other_pisces_tot
#> 1 0 0 0 0 0
#> platichthys_flesus_tot polychaeta_tot saduria_entomon_tot year quarter cruise
#> 1 0 0.05 0 2016 1 BITS
#> predator_spec x y lat long ices_rect pred_size_mm
#> 1 COD 613.9639 6271.507 56.574 16.855 42G6 180
#> pred_weight_g source tot_prey_biom
#> 1 46 Kevin's MSC 20.09
d %>% filter(Unique.pred.id == "2016_1_73_735_COD") %>% as.data.frame()
#> filter: removed 21,729 rows (>99%), 2 rows remaining
#> SubDiv Year Month Day N.with.food N.regurgitated N.skeletal N.empty HN Sample
#> 1 27 2016 2 28 1 NA NA NA 73 735
#> 2 27 2016 2 28 NA NA 1 NA 73 735
#> Length.code Prey.sp.code Prey.size Prey.weight Prey.nr Stage.digestion
#> 1 Lt Clupea harengus 150 19.99 1 1
#> 2 <NA> Bylgides sarsi NA 0.05 NA 2
#> Comment Parasites.in.stomach Quarter Coun. Ship Cruise Pred.size.mm
#> 1 <NA> <NA> 1 SWE NA BITS 180
#> 2 <NA> <NA> 1 SWE NA BITS 180
#> Pred.weight.g Predator.gutted.weight Predator.code Sex Maturity Gall.content
#> 1 46 42 COD M 2 4
#> 2 46 42 COD M 2 4
#> Q.year Age Date Index Lat Long Depth.catch Ices.rect source
#> 1 1_2016 2 <NA> OXBH.2016.73 56.574 16.855 63 42G6 Kevin's MSC
#> 2 1_2016 2 <NA> OXBH.2016.73 56.574 16.855 63 42G6 Kevin's MSC
#> Number Sample.type transect Depthstep Gonad.weight.roundfish
#> 1 NA NA NA NA NA
#> 2 NA NA NA NA NA
#> Liver.weight.roundfish Stomach.content Perc.stomac.content comments Processed
#> 1 NA NA NA NA NA
#> 2 NA NA NA NA NA
#> Unique.pred.id X Y
#> 1 2016_1_73_735_COD 613.9639 6271.507
#> 2 2016_1_73_735_COD 613.9639 6271.507
# Check sample sizes per size class (numbers below bars in Kevins figures)
cod %>%
drop_na(pred_cm_class) %>% # This is because not all predators have length
filter(cruise %in% c("BITS", "BITS-1", "BITS-2")) %>%
filter(year %in% c(2015, 2016, 2017)) %>%
filter(source == "Kevin's MSC") %>%
filter(!pred_cm_class == "(0,6]") %>%
group_by(quarter, pred_cm_class) %>%
summarise(n = n())
#> drop_na: removed 24 rows (<1%), 6,030 rows remaining
#> filter: removed 749 rows (12%), 5,281 rows remaining
#> filter: removed 3,158 rows (60%), 2,123 rows remaining
#> filter: removed 283 rows (13%), 1,840 rows remaining
#> filter: removed 16 rows (1%), 1,824 rows remaining
#> group_by: 2 grouping variables (quarter, pred_cm_class)
#> summarise: now 10 rows and 3 columns, one group variable remaining (quarter)
#> # A tibble: 10 x 3
#> # Groups: quarter [2]
#> quarter pred_cm_class n
#> <dbl> <fct> <int>
#> 1 1 (6,20] 137
#> 2 1 (20,30] 275
#> 3 1 (30,40] 303
#> 4 1 (40,50] 194
#> 5 1 (50,200] 120
#> 6 4 (6,20] 128
#> 7 4 (20,30] 266
#> 8 4 (30,40] 236
#> 9 4 (40,50] 132
#> 10 4 (50,200] 33
# I seem to have more data...
# Check Kevin's raw data...
kev <- read.csv("/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/bentfish/data/stomach_data/Master_Kevin_updated.csv", sep = ";")
head(kev)
#> X coun. ship cruise sample.typ sd ices rect transect date qyear
#> 1 1 SWE DANS BITS trawl SD26-28 43G8 4363 23.11.2015 4_2015
#> 2 2 SWE DANS BITS trawl SD26-28 43G8 4363 23.11.2015 4_2015
#> 3 3 SWE DANS BITS trawl SD26-28 43G8 4363 23.11.2015 4_2015
#> 4 4 SWE DANS BITS trawl SD26-28 43G8 4363 23.11.2015 4_2015
#> 5 5 SWE DANS BITS trawl SD26-28 43G8 4363 23.11.2015 4_2015
#> 6 6 SWE DANS BITS trawl SD26-28 43G8 4363 23.11.2015 4_2015
#> year month day quarter station.haul index lat long depth
#> 1 2015 11 23 4 19 OXBH.2015.19 57,005 18,806 72
#> 2 2015 11 23 4 19 OXBH.2015.19 57,005 18,806 72
#> 3 2015 11 23 4 19 OXBH.2015.19 57,005 18,806 72
#> 4 2015 11 23 4 19 OXBH.2015.19 57,005 18,806 72
#> 5 2015 11 23 4 19 OXBH.2015.19 57,005 18,806 72
#> 6 2015 11 23 4 19 OXBH.2015.19 57,005 18,806 72
#> depthstep sample pred.code pred.size length lengthclass pred.weight
#> 1 >30m 406 COD 190 19 Lc6-20 70
#> 2 >30m 385 COD 320 32 Lc31-40 246
#> 3 >30m 389 COD 330 33 Lc31-40 330
#> 4 >30m 396 COD 300 30 Lc21-30 193
#> 5 >30m 381 COD 310 31 Lc31-40 244
#> 6 >30m 400 COD 250 25 Lc21-30 136
#> pred.gutted.weight age sex maturity gonad.weight liver.weight gall.bladder
#> 1 60 1 M 2 NA 5 1
#> 2 213 3 M 8 NA 9 1
#> 3 499 3 M 2 NA 27 4
#> 4 169 2 M 3 NA 10 4
#> 5 220 2 M 3 NA 11 4
#> 6 115 1 M 3 NA 8 2
#> specimen.note n.full n.regurgitated n.skeletal n.empty stomachcontent
#> 1 1 NA NA NA 0,87
#> 2 1 NA NA NA 1,78
#> 3 NA 1 NA NA 0
#> 4 1 NA NA NA 0,42
#> 5 1 NA NA NA 0,06
#> 6 1 NA NA NA 0,8
#> perc.stomachcontent length.code prey.sp.code prey.size prey.weight prey.nr
#> 1 1,24 Mysis mixta NA 0,87 19
#> 2 0,72 Ls Clupea harengus 85 1,78 1
#> 3 0 NA 0
#> 4 0,22 Mysis mixta NA 0,42 10
#> 5 0,02 Mysis mixta NA 0,06 2
#> 6 0,59 Mysis mixta NA 0,8 19
#> stage.digestion comment parasites.in.stomach processed
#> 1 1 original
#> 2 1 without head original
#> 3 NA fresh scales sure
#> 4 1 original
#> 5 1 original
#> 6 1 original
kev$prey.weight <- as.numeric(gsub(",", ".", gsub("\\.", "", kev$prey.weight)))
kev %>%
filter(lengthclass == "Lc6-20" & pred.code == "COD" & year %in% c(2015, 2016, 2017) & quarter == 1) %>%
drop_na(prey.sp.code) %>%
drop_na(prey.weight) %>%
group_by(prey.sp.code) %>%
summarise(tot_weight_by_prey = sum(prey.weight)) %>%
arrange(desc(tot_weight_by_prey)) %>%
as.data.frame()
#> filter: removed 15,869 rows (99%), 223 rows remaining
#> drop_na: no rows removed
#> drop_na: no rows removed
#> group_by: one grouping variable (prey.sp.code)
#> summarise: now 23 rows and 2 columns, ungrouped
#> prey.sp.code tot_weight_by_prey
#> 1 Clupea harengus 19.99
#> 2 Mysis mixta 6.52
#> 3 Bylgides sarsi 4.52
#> 4 Saduria entomon 3.75
#> 5 Gobiidae 3.53
#> 6 Clupeidae 1.29
#> 7 Diastylis rathkei 1.11
#> 8 Crangon crangon 0.58
#> 9 Priapulida 0.40
#> 10 Gammarus sp. 0.28
#> 11 Mysidae 0.22
#> 12 Halicryptus spinulosus 0.17
#> 13 Cumacea 0.15
#> 14 Limecola balthica 0.11
#> 15 Bivalvia 0.10
#> 16 Pontoporeia femorata 0.07
#> 17 Remains 0.06
#> 18 Amphipoda 0.01
#> 19 Monoporeia affinis 0.01
#> 20 Mytilus sp. 0.01
#> 21 Neomysis integer 0.01
#> 22 0.00
#> 23 Copepoda 0.00
# OK so probably Kevin made some additional filter, because I also have a larger sample sizePlot potential response variables… Unlike Haase, I want to work with individual stomachs
# Plot and calculate proportion empty stomachs
head(data.frame(cod))
#> unique_pred_id amphipoda_tot
#> 1 2007_4_101_2007_11_7_101_DAN2_40G5_1_M_490_1028_1_COD 0.00
#> 2 2007_4_101_2007_11_7_101_DAN2_40G5_118_F_210_82_1_COD 0.24
#> 3 2007_4_101_2007_11_7_101_DAN2_40G5_45_F_480_784_1_COD 0.00
#> 4 2007_4_101_2007_11_7_101_DAN2_40G5_65_F_920_7140_1_COD 0.00
#> 5 2007_4_101_2007_11_7_101_DAN2_40G5_66_M_740_3812_1_COD 0.00
#> 6 2007_4_101_2007_11_7_101_DAN2_40G5_70_F_590_2198_1_COD 0.00
#> bivalvia_tot clupeidae_tot sprattus_sprattus_tot clupea_harengus_tot
#> 1 0 0 5.49 0.00
#> 2 0 0 0.00 0.00
#> 3 0 0 0.00 0.00
#> 4 0 0 26.59 179.66
#> 5 0 0 0.00 88.51
#> 6 0 0 27.32 163.47
#> gadiformes_tot gobiidae_tot mysidae_tot non_bio_tot other_crustacea_tot
#> 1 0 0 0 0 0
#> 2 0 0 0 0 0
#> 3 0 0 0 0 0
#> 4 0 0 0 0 0
#> 5 0 0 0 0 0
#> 6 0 0 0 0 0
#> other_tot other_pisces_tot platichthys_flesus_tot polychaeta_tot
#> 1 0 0 0 0.00
#> 2 0 0 0 0.24
#> 3 0 0 0 0.00
#> 4 0 0 0 0.00
#> 5 0 0 0 0.00
#> 6 0 0 0 0.00
#> saduria_entomon_tot year quarter cruise predator_spec x y
#> 1 0.0 2007 4 BITS-2 COD 520.8793 6176.115
#> 2 0.0 2007 4 BITS-2 COD 520.8793 6176.115
#> 3 0.9 2007 4 BITS-2 COD 520.8793 6176.115
#> 4 0.0 2007 4 BITS-2 COD 520.8793 6176.115
#> 5 0.0 2007 4 BITS-2 COD 520.8793 6176.115
#> 6 0.0 2007 4 BITS-2 COD 520.8793 6176.115
#> lat long ices_rect pred_size_mm pred_weight_g source
#> 1 55.73032 15.33247 40G5 490 1028 Kevin's MSC
#> 2 55.73032 15.33247 40G5 210 82 Kevin's MSC
#> 3 55.73032 15.33247 40G5 480 784 Kevin's MSC
#> 4 55.73032 15.33247 40G5 920 7140 Kevin's MSC
#> 5 55.73032 15.33247 40G5 740 3812 Kevin's MSC
#> 6 55.73032 15.33247 40G5 590 2198 Kevin's MSC
#> tot_prey_biom pred_cm pred_cm_class
#> 1 5.49 49 (40,50]
#> 2 0.48 21 (20,30]
#> 3 0.90 48 (40,50]
#> 4 206.25 92 (50,200]
#> 5 88.51 74 (50,200]
#> 6 190.79 59 (50,200]
cod %>%
mutate(empty_stomach = ifelse(tot_prey_biom == 0, "Y", "N")) %>%
ggplot(., aes(empty_stomach, fill = empty_stomach)) +
geom_bar()
#> mutate: new variable 'empty_stomach' (character) with 2 unique values and 0% NA
# Plot "stomach" fullness, i.e., the weight of prey in stomach (relative to predator weight)
# Can also be called Feeding Ratio, FR
t <- cod %>% drop_na(pred_weight_g)
#> drop_na: removed 1,255 rows (21%), 4,799 rows remaining
t <- cod %>% drop_na(pred_cm)
#> drop_na: removed 24 rows (<1%), 6,030 rows remaining
# Next find the most abundance prey groups for cod and flounder (to see if the density has any effect
# on the biomass of common prey in their stomachs)
# To do that, I need long format again, group by prey item and summarise
colnames(cod)
#> [1] "unique_pred_id" "amphipoda_tot" "bivalvia_tot"
#> [4] "clupeidae_tot" "sprattus_sprattus_tot" "clupea_harengus_tot"
#> [7] "gadiformes_tot" "gobiidae_tot" "mysidae_tot"
#> [10] "non_bio_tot" "other_crustacea_tot" "other_tot"
#> [13] "other_pisces_tot" "platichthys_flesus_tot" "polychaeta_tot"
#> [16] "saduria_entomon_tot" "year" "quarter"
#> [19] "cruise" "predator_spec" "x"
#> [22] "y" "lat" "long"
#> [25] "ices_rect" "pred_size_mm" "pred_weight_g"
#> [28] "source" "tot_prey_biom" "pred_cm"
#> [31] "pred_cm_class"
cod_important_prey <- cod %>%
pivot_longer(2:16) %>%
group_by(name) %>%
summarise(tot_prey = sum(value)) %>%
mutate(percent = round(tot_prey / sum(tot_prey), digits = 2)*100) %>% # calculate also percent of total
arrange(desc(tot_prey)) %>%
mutate(spec = "cod")
#> pivot_longer: reorganized (amphipoda_tot, bivalvia_tot, clupeidae_tot, sprattus_sprattus_tot, clupea_harengus_tot, …) into (name, value) [was 6054x31, now 90810x18]
#> group_by: one grouping variable (name)
#> summarise: now 15 rows and 2 columns, ungrouped
#> mutate: new variable 'percent' (double) with 7 unique values and 0% NA
#> mutate: new variable 'spec' (character) with one unique value and 0% NA
# Same for flounder
fle_important_prey <- fle %>%
pivot_longer(2:16) %>%
group_by(name) %>%
summarise(tot_prey = sum(value)) %>%
mutate(percent = round(tot_prey / sum(tot_prey), digits = 2)*100) %>% # calculate also percent of total
arrange(desc(tot_prey)) %>%
mutate(spec = "fle")
#> pivot_longer: reorganized (amphipoda_tot, bivalvia_tot, clupeidae_tot, sprattus_sprattus_tot, clupea_harengus_tot, …) into (name, value) [was 2257x31, now 33855x18]
#> group_by: one grouping variable (name)
#> summarise: now 15 rows and 2 columns, ungrouped
#> mutate: new variable 'percent' (double) with 9 unique values and 0% NA
#> mutate: new variable 'spec' (character) with one unique value and 0% NA
plotdat <- bind_rows(cod_important_prey, fle_important_prey)
ggplot(plotdat,
aes(name, percent, percent, fill = spec)) +#aes(fct_reorder(name, percent, .desc = TRUE), percent, fill = spec)) +
geom_bar(stat = "identity") +
theme_classic(base_size = 16) +
#coord_flip() +
theme(axis.text.x = element_text(angle = 90)) +
scale_fill_brewer(palette = "Set2") +
NULL
# Ok, so the species that both prey feed on (even though very little(!)) are:
common_prey <- c("amphipoda_tot", "clupea_harengus_tot", "clupeidae_tot", "other_crustacea_tot",
"other_pisces_tot", "polychaeta_tot", "saduria_entomon_tot", "sprattus_sprattus_tot")
# Total feeding ratio, proportion of saduria and proportion of common prey!
# First for cod
cod_fr <- cod %>%
mutate(pred_weight_g = replace_na(pred_weight_g, -9)) %>%
mutate(pred_weight_g = ifelse(pred_weight_g == -9, 0.01*pred_cm^3, pred_weight_g)) %>%
mutate(FR_tot = (tot_prey_biom)/(pred_weight_g - tot_prey_biom),
FR_saduria = (saduria_entomon_tot)/(pred_weight_g - tot_prey_biom)) %>%
# mutate(FR_tot = tot_prey_biom/(pred_weight_g - tot_prey_biom),
# FR_saduria = saduria_entomon_tot/(pred_weight_g - tot_prey_biom)) %>%
mutate(prop_saduria = saduria_entomon_tot/tot_prey_biom,
prop_common = (amphipoda_tot + clupea_harengus_tot + clupeidae_tot + other_crustacea_tot +
other_pisces_tot + polychaeta_tot + saduria_entomon_tot + sprattus_sprattus_tot) / tot_prey_biom,
common_tot = (amphipoda_tot + clupea_harengus_tot + clupeidae_tot + other_crustacea_tot +
other_pisces_tot + polychaeta_tot + saduria_entomon_tot + sprattus_sprattus_tot)) %>%
filter(quarter %in% c(1, 4)) %>%
filter(lat < 58 & lat > 55) %>%
#filter(FR_tot > 0 & FR_tot < 1) #%>% # strange id: 2007_4_72_2007_11_6_72_DAN2_40G4_6_F_250_134_1_COD
filter(FR_tot > 0)
#> mutate: changed 1,255 values (21%) of 'pred_weight_g' (1255 fewer NA)
#> mutate: changed 1,255 values (21%) of 'pred_weight_g' (24 new NA)
#> mutate: new variable 'FR_tot' (double) with 4,821 unique values and <1% NA
#> new variable 'FR_saduria' (double) with 1,229 unique values and <1% NA
#> mutate: new variable 'prop_saduria' (double) with 836 unique values and 12% NA
#> new variable 'prop_common' (double) with 1,242 unique values and 12% NA
#> new variable 'common_tot' (double) with 2,086 unique values and 0% NA
#> filter: removed 91 rows (2%), 5,963 rows remaining
#> filter: removed 20 rows (<1%), 5,943 rows remaining
#> filter: removed 722 rows (12%), 5,221 rows remaining
# I guess in theory the total stomach content could weigh more than the fish without food
# dplyr::select(unique_pred_id, FR_tot) %>%
# arrange(desc(FR_tot)) %>%
# as.data.frame()
# Now flounder
# First for cod
fle_fr <- fle %>%
mutate(pred_weight_g = replace_na(pred_weight_g, -9)) %>%
mutate(pred_weight_g = ifelse(pred_weight_g == -9, 0.01*pred_cm^3, pred_weight_g)) %>%
mutate(FR_tot = (tot_prey_biom)/(pred_weight_g - tot_prey_biom),
FR_saduria = (saduria_entomon_tot)/(pred_weight_g - tot_prey_biom)) %>%
# mutate(FR_tot = tot_prey_biom/(pred_weight_g - tot_prey_biom),
# FR_saduria = saduria_entomon_tot/(pred_weight_g - tot_prey_biom)) %>%
mutate(prop_saduria = saduria_entomon_tot/tot_prey_biom,
prop_common = (amphipoda_tot + clupea_harengus_tot + clupeidae_tot + other_crustacea_tot +
other_pisces_tot + polychaeta_tot + saduria_entomon_tot + sprattus_sprattus_tot) / tot_prey_biom,
common_tot = (amphipoda_tot + clupea_harengus_tot + clupeidae_tot + other_crustacea_tot +
other_pisces_tot + polychaeta_tot + saduria_entomon_tot + sprattus_sprattus_tot)) %>%
filter(quarter %in% c(1, 4)) %>%
filter(lat < 58 & lat > 55) %>%
#filter(FR_tot > 0 & FR_tot < 1) #%>% # strange id: 2007_4_72_2007_11_6_72_DAN2_40G4_6_F_250_134_1_COD
filter(FR_tot > 0)
#> mutate: changed 918 values (41%) of 'pred_weight_g' (918 fewer NA)
#> mutate: changed 918 values (41%) of 'pred_weight_g' (7 new NA)
#> mutate: new variable 'FR_tot' (double) with 1,601 unique values and <1% NA
#> new variable 'FR_saduria' (double) with 640 unique values and <1% NA
#> mutate: new variable 'prop_saduria' (double) with 396 unique values and 25% NA
#> new variable 'prop_common' (double) with 694 unique values and 25% NA
#> new variable 'common_tot' (double) with 495 unique values and 0% NA
#> filter: removed 58 rows (3%), 2,199 rows remaining
#> filter: removed 10 rows (<1%), 2,189 rows remaining
#> filter: removed 522 rows (24%), 1,667 rows remaining
# I guess in theory the total stomach content could weigh more than the fish without food
# dplyr::select(unique_pred_id, FR_tot) %>%
# arrange(desc(FR_tot)) %>%
# as.data.frame()
# What is the mean FR per year
cod %>%
mutate(pred_weight_g = replace_na(pred_weight_g, -9)) %>%
mutate(pred_weight_g = ifelse(pred_weight_g == -9, 0.01*pred_cm^3, pred_weight_g)) %>%
mutate(FR_tot = (100 * tot_prey_biom)/(pred_weight_g - tot_prey_biom),
FR_saduria = (100 * saduria_entomon_tot)/(pred_weight_g - tot_prey_biom)) %>%
filter(FR_tot > -1) %>%
group_by(year) %>%
summarise(mean_FR_tot = mean(FR_tot)) %>%
ggplot(., aes(year, mean_FR_tot)) +
geom_point(size = 4) +
stat_smooth() +
theme_classic(base_size = 18)
#> mutate: changed 1,255 values (21%) of 'pred_weight_g' (1255 fewer NA)
#> mutate: changed 1,255 values (21%) of 'pred_weight_g' (24 new NA)
#> mutate: new variable 'FR_tot' (double) with 4,809 unique values and <1% NA
#> new variable 'FR_saduria' (double) with 1,228 unique values and <1% NA
#> filter: removed 25 rows (<1%), 6,029 rows remaining
#> group_by: one grouping variable (year)
#> summarise: now 13 rows and 2 columns, ungrouped
#> `geom_smooth()` using method = 'loess' and formula 'y ~ x'
# Plot the distribution of FR
cod_fr %>%
ggplot(., aes(FR_tot)) +
geom_histogram() +
theme_classic(base_size = 14) +
coord_cartesian(expand = 0) +
NULL
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
plot_map_raster +
geom_point(data = cod_fr, aes(x = x * 1000, y = y * 1000, color = log(FR_tot))) +
scale_color_viridis() +
facet_wrap(cruise ~ quarter) +
theme_plot()
# Instead plot the proportion of saduria
plot_map_raster +
geom_point(data = cod_fr, aes(x = x * 1000, y = y * 1000, color = prop_saduria)) +
scale_color_viridis() +
#facet_wrap(cruise ~ quarter) +
theme_plot()
# Plot the distribution of proportions - color if the prop is zero!
cod_fr %>%
mutate(saduria_category = ifelse(prop_saduria == 0, "No", "Some")) %>%
mutate(saduria_category = ifelse(prop_saduria == 1, "Only", saduria_category)) %>%
ggplot(., aes(prop_saduria, fill = saduria_category)) +
geom_histogram() +
scale_fill_brewer(palette = "Set2") +
theme_classic(base_size = 14) +
coord_cartesian(expand = 0) +
annotate("text", label = paste("Mean proportion Saduria = ",
round(mean(cod_fr$prop_saduria), digits = 3)), size = 6,
x = 0.5, y = 2000) + # add the mean proportion as text!
NULL
#> mutate: new variable 'saduria_category' (character) with 2 unique values and 0% NA
#> mutate: changed 273 values (5%) of 'saduria_category' (0 new NA)
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
cod_fr %>%
mutate(saduria_category = ifelse(prop_saduria == 0, "No", "Some")) %>%
mutate(saduria_category = ifelse(prop_saduria == 1, "Only", saduria_category)) %>%
ggplot(., aes(prop_saduria, fill = saduria_category)) +
geom_histogram() +
facet_wrap(~ year, scales = "free") +
scale_fill_brewer(palette = "Set2") +
theme_classic(base_size = 14) +
coord_cartesian(expand = 0) +
NULL
#> mutate: new variable 'saduria_category' (character) with 2 unique values and 0% NA
#> mutate: changed 273 values (5%) of 'saduria_category' (0 new NA)
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Conduct test analyses
How does the proportion and total weight of prey per predator weight (prey here being total prey biomass, saduria biomass and biomass of common prey) relate to flounder, cod and total density of the two species? Do the same for flounder!
# First read the density models and predict the values
qwraps2::lazyload_cache_dir(path = "R/analysis/density_spatial_trend_models_cache/html")
#> Lazyloading: R/analysis/density_spatial_trend_models_cache/html/condition model with spatial trend
#> Lazyloading: R/analysis/density_spatial_trend_models_cache/html/fit cod covariates and field
#> Lazyloading: R/analysis/density_spatial_trend_models_cache/html/fit cod depth covariate
#> Lazyloading: R/analysis/density_spatial_trend_models_cache/html/fit cod no covars
#> Lazyloading: R/analysis/density_spatial_trend_models_cache/html/fit fle covariates and field
#> Lazyloading: R/analysis/density_spatial_trend_models_cache/html/fit fle depth covariate
#> Lazyloading: R/analysis/density_spatial_trend_models_cache/html/fit fle no covars
#> Lazyloading: R/analysis/density_spatial_trend_models_cache/html/make barrier spde mesh
#> Lazyloading: R/analysis/density_spatial_trend_models_cache/html/predict on grid
#> Lazyloading: R/analysis/density_spatial_trend_models_cache/html/spatial trends
# Make predictions onto the stomach data defined earlier (1 row per stomach and prey?)
# Add depth and rename coordinates to match the data used for fitting the density model
# Read the tifs
west <- raster("/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/cod_condition/data/depth_geo_tif/D5_2018_rgb-1.tif")
#plot(west)
east <- raster("/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/cod_condition/data/depth_geo_tif/D6_2018_rgb-1.tif")
# plot(east)
dep_rast <- raster::merge(west, east)
cod_fr$depth <- extract(dep_rast, cod_fr[, 24:23])
fle_fr$depth <- extract(dep_rast, fle_fr[, 24:23])
# Convert to depth (instead of elevation)
ggplot(cod_fr, aes(depth)) + geom_histogram()
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
cod_fr$depth <- (cod_fr$depth - max(drop_na(cod_fr)$depth)) *-1
#> drop_na: no rows removed
ggplot(cod_fr, aes(depth)) + geom_histogram()
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
fle_fr$depth <- (fle_fr$depth - max(drop_na(fle_fr)$depth)) *-1
#> drop_na: no rows removed
ggplot(cod_fr, aes(x, y, color = depth)) +
geom_point() +
scale_color_viridis()
cod_fr <- cod_fr %>% rename("X" = "x", "Y" = "y")
#> rename: renamed 2 variables (X, Y)
cod_fr <- cod_fr %>% mutate(year = as.integer(year))
#> mutate: converted 'year' from double to integer (0 new NA)
cod_fr <- cod_fr %>% mutate(year = as.integer(year))
#> mutate: no changes
fle_fr <- fle_fr %>% rename("X" = "x", "Y" = "y")
#> rename: renamed 2 variables (X, Y)
fle_fr <- fle_fr %>% mutate(year = as.integer(year))
#> mutate: converted 'year' from double to integer (0 new NA)
fle_fr <- fle_fr %>% mutate(year = as.integer(year))
#> mutate: no changes
# Standardize depth to the DENSITY data, not this data set
density_dat <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/mdat_cpue.csv")
#>
#> ── Column specification ────────────────────────────────────────────────────────
#> cols(
#> density = col_double(),
#> year = col_double(),
#> lat = col_double(),
#> lon = col_double(),
#> quarter = col_double(),
#> haul.id = col_character(),
#> IDx = col_character(),
#> ices_rect = col_character(),
#> SD = col_double(),
#> density_fle = col_double(),
#> depth = col_double(),
#> oxy = col_double(),
#> temp = col_double(),
#> X = col_double(),
#> Y = col_double()
#> )
cod_fr <- cod_fr %>% mutate(depth_sc = (depth - mean(density_dat$depth)) / sd(density_dat$depth))
#> mutate: new variable 'depth_sc' (double) with 98 unique values and 0% NA
fle_fr <- fle_fr %>% mutate(depth_sc = (depth - mean(density_dat$depth)) / sd(density_dat$depth))
#> mutate: new variable 'depth_sc' (double) with 77 unique values and 0% NA
cod_fr2 <- cod_fr %>% filter(year < 2020) # Remove this after I've refitted the models with 2020!
#> filter: removed 258 rows (5%), 4,963 rows remaining
fle_fr2 <- fle_fr %>% filter(year < 2020) # Remove this after I've refitted the models with 2020!
#> filter: removed 433 rows (26%), 1,234 rows remaining
# Now predict using the density models
cod_stomach_predcod_density <- predict(mcod2, newdata = dplyr::select(cod_fr2, depth_sc, year, X, Y))
#> Warning: Detected a `s()` smoother. Smoothers are penalized in sdmTMB as
#> of version 0.0.19, but used to be unpenalized.
#> You no longer need to specify `k`, specifing `k` now indicates
#> an upper limit on the basis functions, and the degree of wiggliness
#> is determined by the data.
cod_stomach_predfle_density <- predict(mfle2, newdata = dplyr::select(cod_fr2, depth_sc, year, X, Y))
#> Warning: Detected a `s()` smoother. Smoothers are penalized in sdmTMB as
#> of version 0.0.19, but used to be unpenalized.
#> You no longer need to specify `k`, specifing `k` now indicates
#> an upper limit on the basis functions, and the degree of wiggliness
#> is determined by the data.
fle_stomach_predcod_density <- predict(mcod2, newdata = dplyr::select(fle_fr2, depth_sc, year, X, Y))
#> Warning: Detected a `s()` smoother. Smoothers are penalized in sdmTMB as
#> of version 0.0.19, but used to be unpenalized.
#> You no longer need to specify `k`, specifing `k` now indicates
#> an upper limit on the basis functions, and the degree of wiggliness
#> is determined by the data.
fle_stomach_predfle_density <- predict(mfle2, newdata = dplyr::select(fle_fr2, depth_sc, year, X, Y))
#> Warning: Detected a `s()` smoother. Smoothers are penalized in sdmTMB as
#> of version 0.0.19, but used to be unpenalized.
#> You no longer need to specify `k`, specifing `k` now indicates
#> an upper limit on the basis functions, and the degree of wiggliness
#> is determined by the data.
# head(dplyr::select(cod_fr2, X, Y))
# head(predcod_density)
# head(predfle_density)
cod_fr2$predcod_density <- exp(cod_stomach_predcod_density$est)
cod_fr2$predfle_density <- exp(cod_stomach_predfle_density$est)
fle_fr2$predcod_density <- exp(fle_stomach_predcod_density$est)
fle_fr2$predfle_density <- exp(fle_stomach_predfle_density$est)
# Add scales predicted densities to the stomach data
cod_fr2 <- cod_fr2 %>%
mutate(predcod_density_sc = (predcod_density - mean(predcod_density))/sd(predcod_density),
predfle_density_sc = (predfle_density - mean(predfle_density))/sd(predfle_density))
#> mutate: new variable 'predcod_density_sc' (double) with 396 unique values and 0% NA
#> new variable 'predfle_density_sc' (double) with 397 unique values and 0% NA
fle_fr2 <- fle_fr2 %>%
mutate(predcod_density_sc = (predcod_density - mean(predcod_density))/sd(predcod_density),
predfle_density_sc = (predfle_density - mean(predfle_density))/sd(predfle_density))
#> mutate: new variable 'predcod_density_sc' (double) with 168 unique values and 0% NA
#> new variable 'predfle_density_sc' (double) with 168 unique values and 0% NA
# Now plot stomach summaries against predicted densities
head(cod_fr2)
#> # A tibble: 6 x 42
#> unique_pred_id amphipoda_tot bivalvia_tot clupeidae_tot sprattus_spratt…
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 2007_4_101_20… 0 0 0 5.49
#> 2 2007_4_101_20… 0.24 0 0 0
#> 3 2007_4_101_20… 0 0 0 0
#> 4 2007_4_101_20… 0 0 0 26.6
#> 5 2007_4_101_20… 0 0 0 0
#> 6 2007_4_101_20… 0 0 0 27.3
#> # … with 37 more variables: clupea_harengus_tot <dbl>, gadiformes_tot <dbl>,
#> # gobiidae_tot <dbl>, mysidae_tot <dbl>, non_bio_tot <dbl>,
#> # other_crustacea_tot <dbl>, other_tot <dbl>, other_pisces_tot <dbl>,
#> # platichthys_flesus_tot <dbl>, polychaeta_tot <dbl>,
#> # saduria_entomon_tot <dbl>, year <int>, quarter <dbl>, cruise <chr>,
#> # predator_spec <chr>, X <dbl>, Y <dbl>, lat <dbl>, long <dbl>,
#> # ices_rect <chr>, pred_size_mm <dbl>, pred_weight_g <dbl>, source <chr>,
#> # tot_prey_biom <dbl>, pred_cm <dbl>, pred_cm_class <fct>, FR_tot <dbl>,
#> # FR_saduria <dbl>, prop_saduria <dbl>, prop_common <dbl>, common_tot <dbl>,
#> # depth <dbl>, depth_sc <dbl>, predcod_density <dbl>, predfle_density <dbl>,
#> # predcod_density_sc <dbl>, predfle_density_sc <dbl>Plot cod stomach against density
# Common prey
p1 <- ggplot(cod_fr2, aes(predcod_density_sc, prop_common)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2) +
theme_classic(base_size = 13) +
ggtitle("prop common prey vs cod density")
p2 <- ggplot(cod_fr2, aes(predcod_density_sc, common_tot/pred_weight_g)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2) +
theme_classic(base_size = 13) +
ggtitle("g/g common prey vs cod density")
p3 <- ggplot(cod_fr2, aes(predfle_density_sc, prop_common)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2) +
theme_classic(base_size = 13) +
ggtitle("prop common prey vs fle density")
p4 <- ggplot(cod_fr2, aes(predfle_density_sc, common_tot/pred_weight_g)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2) +
theme_classic(base_size = 13) +
ggtitle("g/g common prey vs fle density")
p5 <- ggplot(cod_fr2, aes(predfle_density_sc + predcod_density_sc, prop_common)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2) +
theme_classic(base_size = 13) +
ggtitle("prop common prey vs fle + cod density")
p6 <- ggplot(cod_fr2, aes(predfle_density_sc + predcod_density_sc, common_tot/pred_weight_g)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2) +
theme_classic(base_size = 13) +
ggtitle("g/g common prey vs fle + coddensity")
(p1|p2) / (p3|p4) / (p5|p6) + plot_annotation(title = "Cod stomachs",
theme = theme(plot.title = element_text(size = 16)))
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
# Saduria
p7 <- ggplot(cod_fr2, aes(predcod_density_sc, prop_saduria)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2) +
theme_classic(base_size = 13) +
ggtitle("prop saduria vs cod density")
p8 <- ggplot(cod_fr2, aes(predcod_density_sc, saduria_entomon_tot/pred_weight_g)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2) +
theme_classic(base_size = 13) +
ggtitle("g/g saduria vs cod density")
p9 <- ggplot(cod_fr2, aes(predfle_density_sc, prop_saduria)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2) +
theme_classic(base_size = 13) +
ggtitle("prop saduria vs fle density")
p10 <- ggplot(cod_fr2, aes(predfle_density_sc, saduria_entomon_tot/pred_weight_g)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2) +
theme_classic(base_size = 13) +
ggtitle("g/g saduria vs fle density")
p11 <- ggplot(cod_fr2, aes(predfle_density_sc + predcod_density_sc, prop_saduria)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2) +
theme_classic(base_size = 13) +
ggtitle("prop saduria vs fle + cod density")
p12 <- ggplot(cod_fr2, aes(predfle_density_sc + predcod_density_sc, saduria_entomon_tot/pred_weight_g)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2) +
theme_classic(base_size = 13) +
ggtitle("g/g saduria vs fle + coddensity")
(p7|p8) / (p9|p10) / (p11|p12) + plot_annotation(title = "Cod stomachs",
theme = theme(plot.title = element_text(size = 16)))
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
# Total stomach content (FR)
p13 <- ggplot(cod_fr2, aes(predcod_density_sc, FR_tot)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2) +
theme_classic(base_size = 13) +
ggtitle("FR vs cod density")
p14 <- ggplot(cod_fr2, aes(predfle_density_sc, FR_tot)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2) +
theme_classic(base_size = 13) +
ggtitle("FR vs fle density")
p15 <- ggplot(cod_fr2, aes(predfle_density_sc + predcod_density, FR_tot)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2) +
theme_classic(base_size = 13) +
ggtitle("FR vs fle + coddensity")
(p13/p14/p15) + plot_annotation(title = "Cod stomachs", theme = theme(plot.title = element_text(size = 16)))
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
# Closer inspection of some of the plots...
ggplot(cod_fr2, aes(predfle_density_sc, prop_saduria)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2, method = "lm") +
theme_classic(base_size = 13) +
ggtitle("prop saduria vs fle density")
#> `geom_smooth()` using formula 'y ~ x'
summary(lm(cod_fr2$prop_saduria ~ cod_fr2$predfle_density))
#>
#> Call:
#> lm(formula = cod_fr2$prop_saduria ~ cod_fr2$predfle_density)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -0.1330 -0.1321 -0.1271 -0.1001 0.9968
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 1.330e-01 4.423e-03 30.082 < 2e-16 ***
#> cod_fr2$predfle_density -1.882e-04 3.076e-05 -6.119 1.02e-09 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 0.2942 on 4961 degrees of freedom
#> Multiple R-squared: 0.00749, Adjusted R-squared: 0.00729
#> F-statistic: 37.44 on 1 and 4961 DF, p-value: 1.015e-09
ggplot(cod_fr2, aes(predfle_density, tot_prey_biom/pred_weight_g)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2) +
theme_classic(base_size = 13) +
ggtitle("g/g tot. stomach vs fle density")
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
ggplot(cod_fr2, aes(predfle_density + predcod_density, tot_prey_biom/pred_weight_g)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2) +
theme_classic(base_size = 13) +
ggtitle("g/g tot. stomach vs fle + coddensity")
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
Plot flounder stomach against density
# Common prey
p1 <- ggplot(fle_fr2, aes(predcod_density_sc, prop_common)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2) +
theme_classic(base_size = 13) +
ggtitle("prop common prey vs cod density")
p2 <- ggplot(fle_fr2, aes(predcod_density_sc, common_tot/pred_weight_g)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2) +
theme_classic(base_size = 13) +
ggtitle("g/g common prey vs cod density")
p3 <- ggplot(fle_fr2, aes(predfle_density_sc, prop_common)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2) +
theme_classic(base_size = 13) +
ggtitle("prop common prey vs fle density")
p4 <- ggplot(fle_fr2, aes(predfle_density_sc, common_tot/pred_weight_g)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2) +
theme_classic(base_size = 13) +
ggtitle("g/g common prey vs fle density")
p5 <- ggplot(fle_fr2, aes(predfle_density_sc + predcod_density_sc, prop_common)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2) +
theme_classic(base_size = 13) +
ggtitle("prop common prey vs fle + cod density")
p6 <- ggplot(fle_fr2, aes(predfle_density_sc + predcod_density_sc, common_tot/pred_weight_g)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2) +
theme_classic(base_size = 13) +
ggtitle("g/g common prey vs fle + coddensity")
(p1|p2) / (p3|p4) / (p5|p6) + plot_annotation(title = "Flounder stomachs",
theme = theme(plot.title = element_text(size = 16)))
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
# Saduria
p7 <- ggplot(fle_fr2, aes(predcod_density_sc, prop_saduria)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2) +
theme_classic(base_size = 13) +
ggtitle("prop saduria vs cod density")
p8 <- ggplot(fle_fr2, aes(predcod_density_sc, saduria_entomon_tot/pred_weight_g)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2) +
theme_classic(base_size = 13) +
ggtitle("g/g saduria vs cod density")
p9 <- ggplot(fle_fr2, aes(predfle_density_sc, prop_saduria)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2) +
theme_classic(base_size = 13) +
ggtitle("prop saduria vs fle density")
p10 <- ggplot(fle_fr2, aes(predfle_density_sc, saduria_entomon_tot/pred_weight_g)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2) +
theme_classic(base_size = 13) +
ggtitle("g/g saduria vs fle density")
p11 <- ggplot(fle_fr2, aes(predfle_density_sc + predcod_density_sc, prop_saduria)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2) +
theme_classic(base_size = 13) +
ggtitle("prop saduria vs fle + cod density")
p12 <- ggplot(fle_fr2, aes(predfle_density_sc + predcod_density_sc, saduria_entomon_tot/pred_weight_g)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2) +
theme_classic(base_size = 13) +
ggtitle("g/g saduria vs fle + coddensity")
(p7|p8) / (p9|p10) / (p11|p12) + plot_annotation(title = "Flounder stomachs",
theme = theme(plot.title = element_text(size = 16)))
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
# FR
p13 <- ggplot(fle_fr2, aes(predcod_density_sc, FR_tot)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2) +
theme_classic(base_size = 13) +
ggtitle("FR vs cod density")
p14 <- ggplot(fle_fr2, aes(predfle_density_sc, FR_tot)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2) +
theme_classic(base_size = 13) +
ggtitle("FR vs fle density")
p15 <- ggplot(fle_fr2, aes(predfle_density_sc + predcod_density, FR_tot)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2) +
theme_classic(base_size = 13) +
ggtitle("FR vs fle + coddensity")
(p13/p14/p15) + plot_annotation(title = "Flounder stomachs", theme = theme(plot.title = element_text(size = 16)))
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
More in depth plots…
# Does it seem like proportions of the diet changes with density, but not the total amount of food?
# Check saduria and flounder for instance.
pp7 <- ggplot(fle_fr2, aes(predcod_density_sc, prop_saduria)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2, method = "lm") +
theme_classic(base_size = 13) +
ggtitle("prop saduria vs cod density")
pp8 <- ggplot(filter(fle_fr2, (saduria_entomon_tot/pred_weight_g) < 0.1),
aes(predcod_density_sc, saduria_entomon_tot/pred_weight_g)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2, method = "lm") +
theme_classic(base_size = 13) +
ggtitle("g/g saduria vs cod density")
#> filter: removed 2 rows (<1%), 1,232 rows remaining
pp13 <- ggplot(filter(fle_fr2, (saduria_entomon_tot/pred_weight_g) < 0.1),
aes(predcod_density_sc, FR_tot)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2, method = "lm") +
theme_classic(base_size = 13) +
ggtitle("g/g tot. stomach vs cod density")
#> filter: removed 2 rows (<1%), 1,232 rows remaining
ppfle <- pp7/pp8/pp13 +
plot_layout(ncol = 3) +
plot_annotation(title = "Flounder stomachs", theme = theme(plot.title = element_text(size = 16)))
# And that cod density is unaffected (or less at least)
p7 <- ggplot(filter(cod_fr2, (saduria_entomon_tot/pred_weight_g) < 0.02 & predfle_density_sc < 2),
aes(predfle_density_sc, prop_saduria)) + # Crop it up a little!
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2, method = "lm") +
theme_classic(base_size = 13) +
ggtitle("prop saduria vs fle density")
#> filter: removed 171 rows (3%), 4,792 rows remaining
p8 <- ggplot(filter(cod_fr2, (saduria_entomon_tot/pred_weight_g) < 0.02 & predfle_density_sc < 2),
aes(predfle_density_sc, saduria_entomon_tot/pred_weight_g)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2, method = "lm") +
theme_classic(base_size = 13) +
ggtitle("g/g saduria vs fle density")
#> filter: removed 171 rows (3%), 4,792 rows remaining
p13 <- ggplot(filter(cod_fr2, (saduria_entomon_tot/pred_weight_g) < 0.02 & predfle_density_sc < 2 & FR_tot < 0.5),
aes(predfle_density_sc, FR_tot)) +
geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) +
stat_smooth(size = 2, method = "lm") +
theme_classic(base_size = 13) +
ggtitle("g/g tot. stomach vs fle density")
#> filter: removed 179 rows (4%), 4,784 rows remaining
ppcod <- p7/p8/p13 +
plot_layout(ncol = 3) +
plot_annotation(title = "Cod stomachs", theme = theme(plot.title = element_text(size = 16)))
ppfle / ppcod
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
# What does it mean? Yes, the proportion *may* be negatively affected by densities,
# but the food content is not, which I guess would matter for conditionTry and fit sdmTMB models with densities as covariates…
# Cod
pcod_spde <- make_mesh(cod_fr2, c("X", "Y"), n_knots = 150, type = "kmeans", seed = 42)
plot(pcod_spde)
pfle_spde <- make_mesh(fle_fr2, c("X", "Y"), n_knots = 70, type = "kmeans", seed = 42)
plot(pfle_spde)
# (add small values if 0, subtract small values if 1) for Beta regressions
cod_fr2 <- cod_fr2 %>%
mutate(saduria_entomon_per_mass = saduria_entomon_tot/pred_weight_g,
prop_saduria2 = ifelse(prop_saduria == 0, prop_saduria + 0.0001, prop_saduria),
prop_saduria2 = ifelse(prop_saduria == 1, prop_saduria - 0.0001, prop_saduria2))
#> mutate: new variable 'saduria_entomon_per_mass' (double) with 1,066 unique values and 0% NA
#> new variable 'prop_saduria2' (double) with 782 unique values and 0% NA
fle_fr2 <- fle_fr2 %>%
mutate(saduria_entomon_per_mass = saduria_entomon_tot/pred_weight_g,
prop_saduria2 = ifelse(prop_saduria == 0, prop_saduria + 0.0001, prop_saduria),
prop_saduria2 = ifelse(prop_saduria == 1, prop_saduria - 0.0001, prop_saduria2))
#> mutate: new variable 'saduria_entomon_per_mass' (double) with 430 unique values and 0% NA
#> new variable 'prop_saduria2' (double) with 300 unique values and 0% NA
# OK, proportion models look absolutely AWFUL. I will model stomach content with a tweedie instead
# Proportion Saduria
# mcodsad <- sdmTMB(
# data = cod_fr2,
# formula = prop_saduria2 ~ 0 + as.factor(year) + predfle_density_sc + predcod_density_sc,
# time = "year", fields = "IID", spatial_only = FALSE, spde = pcod_spde, family = student(df = 5))
mcodsad <- sdmTMB(
data = cod_fr2,
formula = saduria_entomon_per_mass ~ 0 + as.factor(year) + s(predfle_density_sc) + s(predcod_density_sc),
time = "year", fields = "IID", spatial_only = FALSE, spde = pcod_spde, family = tweedie())
#> Warning: Detected a `s()` smoother. Smoothers are penalized in sdmTMB as
#> of version 0.0.19, but used to be unpenalized.
#> You no longer need to specify `k`, specifing `k` now indicates
#> an upper limit on the basis functions, and the degree of wiggliness
#> is determined by the data.
cod_fr2$resids_mcodsad <- residuals(mcodsad) # randomized quantile residuals
#> Warning: Detected a `s()` smoother. Smoothers are penalized in sdmTMB as
#> of version 0.0.19, but used to be unpenalized.
#> You no longer need to specify `k`, specifing `k` now indicates
#> an upper limit on the basis functions, and the degree of wiggliness
#> is determined by the data.
hist(cod_fr2$resids_mcodsad)
qqnorm(cod_fr2$resids_mcodsad); abline(a = 0, b = 1)
tidy(mcodsad)
#> term estimate std.error
#> 1 as.factor(year)2007 -7.559057 0.7891209
#> 2 as.factor(year)2008 -7.800466 0.8602433
#> 3 as.factor(year)2009 -8.748290 0.8173057
#> 4 as.factor(year)2010 -9.210658 1.0224320
#> 5 as.factor(year)2011 -10.058870 0.9114980
#> 6 as.factor(year)2012 -8.706427 0.8372932
#> 7 as.factor(year)2013 -9.184026 0.7839526
#> 8 as.factor(year)2015 -8.567464 0.7569644
#> 9 as.factor(year)2016 -9.869242 0.6801726
#> 10 as.factor(year)2017 -10.710367 0.7368018
#> 11 as.factor(year)2018 -9.698110 0.7087329
#> 12 as.factor(year)2019 -10.655589 0.8360574
#> 13 s(predfle_density_sc).1 NA NA
#> 14 s(predfle_density_sc).2 NA NA
#> 15 s(predfle_density_sc).3 NA NA
#> 16 s(predfle_density_sc).4 NA NA
#> 17 s(predfle_density_sc).5 NA NA
#> 18 s(predfle_density_sc).6 NA NA
#> 19 s(predfle_density_sc).7 NA NA
#> 20 s(predfle_density_sc).8 NA NA
#> 21 s(predfle_density_sc).9 NA NA
#> 22 s(predcod_density_sc).1 NA NA
#> 23 s(predcod_density_sc).2 NA NA
#> 24 s(predcod_density_sc).3 NA NA
#> 25 s(predcod_density_sc).4 NA NA
#> 26 s(predcod_density_sc).5 NA NA
#> 27 s(predcod_density_sc).6 NA NA
#> 28 s(predcod_density_sc).7 NA NA
#> 29 s(predcod_density_sc).8 NA NA
#> 30 s(predcod_density_sc).9 NA NA
summary(mcodsad)
#> Spatial model fit by ML ['sdmTMB']
#> Formula: saduria_entomon_per_mass ~ 0 + as.factor(year) + s(predfle_density_sc) +
#> Formula: s(predcod_density_sc)
#> Time column: "year"
#> SPDE: pcod_spde
#> Data: cod_fr2
#> Family: tweedie(link = 'log')
#> coef.est coef.se
#> as.factor(year)2007 -7.56 0.79
#> as.factor(year)2008 -7.80 0.86
#> as.factor(year)2009 -8.75 0.82
#> as.factor(year)2010 -9.21 1.02
#> as.factor(year)2011 -10.06 0.91
#> as.factor(year)2012 -8.71 0.84
#> as.factor(year)2013 -9.18 0.78
#> as.factor(year)2015 -8.57 0.76
#> as.factor(year)2016 -9.87 0.68
#> as.factor(year)2017 -10.71 0.74
#> as.factor(year)2018 -9.70 0.71
#> as.factor(year)2019 -10.66 0.84
#>
#> Matern range parameter: 42.38
#> Dispersion parameter: 0.17
#> Spatial SD: 2.29
#> Spatiotemporal SD: 1.43
#> ML criterion at convergence: -3251.017
#>
#> See ?tidy.sdmTMB to extract these values as a data frame.
mflesad <- sdmTMB(
data = fle_fr2,
formula = saduria_entomon_per_mass ~ 0 + as.factor(year) + s(predfle_density_sc) + s(predcod_density_sc),
time = "year", fields = "IID", spatial_only = FALSE, spde = pfle_spde, family = tweedie())
#> Warning: Detected a `s()` smoother. Smoothers are penalized in sdmTMB as
#> of version 0.0.19, but used to be unpenalized.
#> You no longer need to specify `k`, specifing `k` now indicates
#> an upper limit on the basis functions, and the degree of wiggliness
#> is determined by the data.
tidy(mflesad)
#> term estimate std.error
#> 1 as.factor(year)2015 -7.324439 1.072835
#> 2 as.factor(year)2016 -7.231370 1.020008
#> 3 as.factor(year)2017 -7.521263 1.017771
#> 4 as.factor(year)2018 -7.248929 1.060067
#> 5 as.factor(year)2019 -7.649248 1.130401
#> 6 s(predfle_density_sc).1 NA NA
#> 7 s(predfle_density_sc).2 NA NA
#> 8 s(predfle_density_sc).3 NA NA
#> 9 s(predfle_density_sc).4 NA NA
#> 10 s(predfle_density_sc).5 NA NA
#> 11 s(predfle_density_sc).6 NA NA
#> 12 s(predfle_density_sc).7 NA NA
#> 13 s(predfle_density_sc).8 NA NA
#> 14 s(predfle_density_sc).9 NA NA
#> 15 s(predcod_density_sc).1 NA NA
#> 16 s(predcod_density_sc).2 NA NA
#> 17 s(predcod_density_sc).3 NA NA
#> 18 s(predcod_density_sc).4 NA NA
#> 19 s(predcod_density_sc).5 NA NA
#> 20 s(predcod_density_sc).6 NA NA
#> 21 s(predcod_density_sc).7 NA NA
#> 22 s(predcod_density_sc).8 NA NA
#> 23 s(predcod_density_sc).9 NA NA
summary(mflesad)
#> Spatial model fit by ML ['sdmTMB']
#> Formula: saduria_entomon_per_mass ~ 0 + as.factor(year) + s(predfle_density_sc) +
#> Formula: s(predcod_density_sc)
#> Time column: "year"
#> SPDE: pfle_spde
#> Data: fle_fr2
#> Family: tweedie(link = 'log')
#> coef.est coef.se
#> as.factor(year)2015 -7.32 1.07
#> as.factor(year)2016 -7.23 1.02
#> as.factor(year)2017 -7.52 1.02
#> as.factor(year)2018 -7.25 1.06
#> as.factor(year)2019 -7.65 1.13
#>
#> Matern range parameter: 79.85
#> Dispersion parameter: 0.22
#> Spatial SD: 2.73
#> Spatiotemporal SD: 0.68
#> ML criterion at convergence: -1155.956
#>
#> See ?tidy.sdmTMB to extract these values as a data frame.
fle_fr2$resids_mflesad <- residuals(mflesad) # randomized quantile residuals
#> Warning: Detected a `s()` smoother. Smoothers are penalized in sdmTMB as
#> of version 0.0.19, but used to be unpenalized.
#> You no longer need to specify `k`, specifing `k` now indicates
#> an upper limit on the basis functions, and the degree of wiggliness
#> is determined by the data.
hist(fle_fr2$resids_mflesad)
qqnorm(fle_fr2$resids_mflesad); abline(a = 0, b = 1)
# Plot marginal effects
## Flounder
# Flounder density
nd_fle_fle <- data.frame(predfle_density_sc = seq(min(fle_fr2$predfle_density_sc),
max(fle_fr2$predfle_density_sc), length.out = 100))
nd_fle_fle$year <- 2018L
nd_fle_fle$predcod_density_sc <- 0
pred_fle_fle <- predict(mflesad, newdata = nd_fle_fle, se_fit = TRUE, re_form = NA)
#> Warning: Detected a `s()` smoother. Smoothers are penalized in sdmTMB as
#> of version 0.0.19, but used to be unpenalized.
#> You no longer need to specify `k`, specifing `k` now indicates
#> an upper limit on the basis functions, and the degree of wiggliness
#> is determined by the data.
p1 <- ggplot(pred_fle_fle, aes(predfle_density_sc, exp(est),
ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
geom_line() + geom_ribbon(alpha = 0.4) +
theme_classic(base_size = 14) +
ggtitle("Fle stomach, flounder margin")
# Cod density
nd_fle_cod <- data.frame(predcod_density_sc = seq(min(fle_fr2$predcod_density_sc),
max(fle_fr2$predcod_density_sc), length.out = 100))
nd_fle_cod$year <- 2018L
nd_fle_cod$predfle_density_sc <- 0
pred_fle_cod <- predict(mflesad, newdata = nd_fle_cod, se_fit = TRUE, re_form = NA)
#> Warning: Detected a `s()` smoother. Smoothers are penalized in sdmTMB as
#> of version 0.0.19, but used to be unpenalized.
#> You no longer need to specify `k`, specifing `k` now indicates
#> an upper limit on the basis functions, and the degree of wiggliness
#> is determined by the data.
p2 <- ggplot(pred_fle_cod, aes(predcod_density_sc, exp(est),
ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
geom_line() + geom_ribbon(alpha = 0.4) +
theme_classic(base_size = 14) +
ggtitle("Fle stomach, cod margin")
flemargin <- p1 + p2 + plot_annotation(title = "Flounder marginal effects",
theme = theme(plot.title = element_text(size = 16)))
## Cod
# Flounder density
nd_cod_fle <- data.frame(predfle_density_sc = seq(min(cod_fr2$predfle_density_sc),
max(cod_fr2$predfle_density_sc), length.out = 100))
nd_cod_fle$year <- 2018L
nd_cod_fle$predcod_density_sc <- 0
pred_cod_fle <- predict(mcodsad, newdata = nd_cod_fle, se_fit = TRUE, re_form = NA)
#> Warning: Detected a `s()` smoother. Smoothers are penalized in sdmTMB as
#> of version 0.0.19, but used to be unpenalized.
#> You no longer need to specify `k`, specifing `k` now indicates
#> an upper limit on the basis functions, and the degree of wiggliness
#> is determined by the data.
p3 <- ggplot(pred_cod_fle, aes(predfle_density_sc, exp(est),
ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
geom_line() + geom_ribbon(alpha = 0.4) +
theme_classic(base_size = 14) +
ggtitle("Cod stomach, flounder margin")
# Cod density
nd_cod_cod <- data.frame(predcod_density_sc = seq(min(cod_fr2$predcod_density_sc),
max(cod_fr2$predcod_density_sc), length.out = 100))
nd_cod_cod$year <- 2018L
nd_cod_cod$predfle_density_sc <- 0
pred_cod_cod <- predict(mcodsad, newdata = nd_cod_cod, se_fit = TRUE, re_form = NA)
#> Warning: Detected a `s()` smoother. Smoothers are penalized in sdmTMB as
#> of version 0.0.19, but used to be unpenalized.
#> You no longer need to specify `k`, specifing `k` now indicates
#> an upper limit on the basis functions, and the degree of wiggliness
#> is determined by the data.
p4 <- ggplot(pred_cod_cod, aes(predcod_density_sc, exp(est),
ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
geom_line() + geom_ribbon(alpha = 0.4) +
theme_classic(base_size = 14) +
ggtitle("Cod stomach, cod margin")
(p1 | p2) / (p3 | p4 )